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Abstract 

We study the first-order relativistic correction to the associated production of J /ip with light 
hadrons at B factory experiments at ^/s = 10.58 GeV, in the context of NRQCD factorization. 
We employ a strategy for NRQCD expansion that slightly deviates from the orthodox doctrine, 
in that the matching coefficients are not truly of "short-distance" nature, but explicitly depend 
upon physical kinematic variables rather than partonic ones. Our matching method, with validity 
guaranteed by the Gremm-Kapustin relation, is particularly suited for the inclusive quarkonium 
production and decay processes with involved kinematics, exemplified by the process e^e~ — >■ 
J/il) + gg considered in this work. Despite some intrinsic ambiguity affiliated with the order- 
NRQCD matrix element, if we choose its value as what has been extracted from a recent 
Cornell-potential-model-based analysis, including the relative order-u^ effect is found to increase 
the lowest-order prediction for the integrated J /tp cross section by about 30%, and exert a modest 
impact on J /tjj energy, angular and polarization distributions except near the very upper end of the 
J /il) energy. The order-u'^ contribution to the energy spectrum becomes logarithmically divergent 
at the maximum of J/ip energy. A consistent analysis may require that these large end-point 
logarithms be resummed to all orders in a^. 
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I. INTRODUCTION 



Inclusive production of heavy quarkonium (especially J/^) in various high-energy col- 
lider experiments has long been an intriguing and interesting topic, to which a vast number 
of works have been devoted [1, 2]. To date the mainstream of theoretical investigations 
in this subject is based upon the nonrelativistic QCD (NRQCD) factorization approach, a 
formalism combining the effective-field-theory machinery together with the hard-scattering 
factorization [3] . In the context of NRQCD factorization, the inclusive quarkonium produc- 
tion rate can be expressed in a factorized form, that is, an infinite sum of products of the 
perturbatively calculable partonic cross sections and nonperturbative but universal NRQCD 
matrix elements. One great virtue of this approach is that its predictions can in principle be 
systematically improved. This approach systemizes, and, extends, the conventional color- 
singlet model (CSM). One striking, and, probably also disputable, ingredient of NRQCD 
factorization is the so-called color-octet mechanism, that a heavy quark-antiquark pair in 
a color-octet configuration created in a hard process, is presumed to have non-negligible 
probability to transition into a physical quarkonium state plus additional soft light hadrons. 
Historically, the rapid popularity gained by this novel mechanism is perhaps due to its 
economic explanation of the so-called '■0' surplus puzzle' [4] . 

Although NRQCD factorization has enjoyed considerable successes in many inclusive 
quarkonium decay and production processes, it also faces some serious challenges. Most 
notably, the recent Fermilab Tevatron measurement for J/ ip polarization at large px seems to 
contradict with the benchmark predictions of the color-octet mechanism, i.e. the increasingly 
transverse polarization of the hadro-produced J /il) with increasing pt [5]. Moreover, there is 
also problem for J/ip production in e+e~ collision experiments. The color-octet mechanism 
also anticipates that an enhanced number of J /ip populate near the maximum energy region 
in e'^e~ annihilation, but unfortunately, this quite distinct signature has also not been 
confirmed by recent B factory experiments. 

These acute discrepancies have triggered a great wave of theoretical efforts in recent 
years. Recent technical advancement makes it possible, for the first time, to compute the 
rather involved next-to- leading QCD corrections to J/ip hadroproduction in color-singlet 
channel [6-10], and its effects seem to be quite significant, i.e., to enhance the leading-order 
CSM contribution enormously. This may indicate that, the phenomenological impetus to 
including color-octet contribution seems not as indispensable as that in a decade ago, and 
the correct magnitudes of color-octet matrix elements might be considerably smaller than 
the old numbers extracted by implementing the LO CSM analysis only. 

It is worth emphasizing that, the NRQCD factorization theorems for quarkonium pro- 
duction are only at a conjectural level, which have never been proven rigourously to hold to 
all orders in as- Notwithstanding the great utility of the improvement on the short-distance 
coefficients, it is perhaps more urgent, from the theoretical perspective, to reexamine every 
assumption underlying the nonperturbative aspects of NRQCD factorization, especially for 
the color-octet mechanism. As one of the important progresses along this line, the valid- 
ity of the factorization theorems at two-loop level for gluon-to-quarkonium fragmentation 
function has recently been established after some suitable refinement of the original color- 
octet NRQCD production operator [11] ^. There is also suspect about the applicabifity of 



Now it becomes clear that in some case NRQCD factorization certainly will fail. For example, a novel 
phenomenon dubbed color transfer mechanism [12], was discovered in the production of J/V' comoving 
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the NRQCD velocity-scaling rule to charmoniTim, in particular it was suggested that the 
spin-flip matrix element may play an important role for the hardronization of color-octet cc 
pair [13, 14]. A more serious problem is that since each NRQCD matrix element is a number 
instead of a distribution, so NRQCD factorization makes rather restrictive predictions to 
the various J/'^ energy spectra. It has long been suggested that in certain kinematic region 
of quarkonium production, resummation of a class of enhanced nonperturbative effects is 
crucial to make reliable prediction, which effectively promotes the local NRQCD matrix 
element to a nonperturbative shape function [15, 16]. It is also worth noting that, there is 
also an ongoing endeavor to circumvent the velocity expansion framework of NRQCD, by 
introducing a more general set of fragmentation functions in conventional perturbative QCD 
(pQCD) factorization base to describe inclusive quarkonium production at large pt [17]. 

In recent years, B factories also prove to be another active field for the study of char- 
monium production. The simplicity of the initial e+e~ state, together with the enormous 
integrated luminosity, make the theoretical analysis of the charmonium production process 
particularly clean and fertile. Some recent measurements at B factories have also posed 
challenges to our understanding of charmonium production. One is the unexpectedly large 
cross sections for several exclusive double charmonium production processes in continuum 
e'^e~ annihilation. For example, the cross section for producing J/ip + r)c was first measured 
by the BELLE collaboration [18], which turns out to be almost one order-of-magnitude larger 
than the leading order (LO) NRQCD predictions [19]. After various theoretical works from 
different angles, the consensus now is that after including the large QCD perturbative cor- 
rection [20, 21], in combination with relativistic corrections, this disquieting discrepancy was 
claimed to be largely resolved within the context of NRQCD factorization [22, 23]. 

Another more perplexing observation arises from the inclusive production of J/iJj. The 
production of J/ ip in association with extra charms, is found, quite counter- intuitively, to oc- 
cur much more copiously than that in association with a non-charm final state. The fraction 
of number of events for J/ip plus charmed hadrons to that of the inclusive J/ip events, con- 
ventionally denoted Rcc, was first measured by Belle collaboration to be 0.59lg ;[3±0.12 [18], 
later even shifted to 0.82 ±0.15 ±0.14 [24]. This experimental values are in stark contrast to 
the leading-order (LO) NRQCD predictions to this ratio, which is only about 0.1 ^. Other 
theoretical approaches, e.g., the estimate based on quark-hadron duality hypothesis [32] and 
the color-evaporation model [33] also predict a quite small Rcc- 

Once upon a time, the total J/ip production rate measured at B factories appeared to 
be quite large, i.e., 2.5 pb measured by Babar [34] and 1.5 pb [35] by Belle, which seems 
to request a sizable color-octet contribution such as e+e^ — )■ cc(^Sl^\ "^Pj^^) + g. The color- 
octet effect for J/ip production in e+e" annihilation was first investigated by Braaten and 
Chen [36] (see also [37], and for a very recent study of the NLO perturbative correction, see 
[38].). However, including this contribution will further dilute the ratio Rcc- Furthermore, 
an unusual signature of this mechanism is that an end-point peak is expected in the J/tjj 
energy spectrum. Unfortunately, there is no experimental evidences for the existence of 
such a peak [34, 35]. To rescue the color-octet mechanism, later on the end-point Sudakov 



with an additional heavy quark. In this case, soft color exchange between the comoving quark and the 
constitutes of J/tl) may invalidate the NRQCD factorization at two-loop level. 
^ Note that the LO NRQCD predictions to the Jf-ip production associated with charmed or noncharmed 
hadrons are identical to the CSM predictions [25-31], i.e. to proceed through the parton processes 
e+e~ J/ip + CC and e+e" J/ip + gg, respectively. 
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logarithms have been identified and resummed, together with introduction of a phenomeno- 
logical shape function, one can show that the J/ip energy distribution can be smeared out 
in accordance with the data [39] . It is worth mentioning that absolute normalization of the 
color-octet contribution is not affected much by including these refinements, and by that 
time its contribution was assumed to predominate over the color-singlet contribution. 

In the past couple of years, significant progresses toward resolving these puzzles have 
been made from both experimental and theoretical angles. From the experimental side, 
recently Belle collaboration was able to precisely measure the cross sections for prompt 
J/%1) production in association with charmed and non-charmed states separately [40]: 



This new measurement has not subtracted the feeddown contribution from il){2S). 

The most important recent theoretical progress in this subject is perhaps the fulfillment of 
the NLO QCD corrections to both channels. It turns out that the inclusion of the NLO QCD 
correction significantly enhances the J/il) + cc production rate [41, 42]. Recently, the next-to- 
leading (NLO) perturbative correction to e+e~ — )■ J/ip + gg has also been conducted, which 
enhances the LO cross section by only about 20% [43, 44] ^ . The significant enhancement 
to the former and the modest one to the latter is of help for the predicted i?cc value to 
approach the measured one. When the feeddown effects are included, the rough agreement 
seems also to be achieved for both the associated J/'4> production subprocesses, so the 
alarming discrepancies seem to be greatly alleviated ^. 

The above analysis tends to indicate that, CSM alone seems sufficient to explain the 
data, and there seems no much room left for the color-octet contribution. As a result, the 
color-octet matrix element may be considerably smaller than what was used to be assumed 
when fitted from the Tevatron data. Nevertheless, it is still premature to assert that we 
already have satisfactory understanding of inclusive J/ip production at B factory, because 
there is still one important component missing. That is, in compliance with the NRQCD 
power counting, one should also take the first-order relativistic correction into account, 
since its effect is parametrically more important than the color-octet contribution. This 
is particularly relevant for e^e~ J /if) + Xng^t since the size of NLO QCD correction 
to e^e~ — > J/ip + gg is mild ^. It is thus interesting to examine whether the relativistic 
correction brings in sizable effects to this process or not. 

The main purpose of this work is to answer this question, that is, to calculate the first- 
order relativistic correction to the inclusive J/ip production rate in e'^e" J/ijj + -'^light in 
NRQCD factorization. Concretely, we will be considering the process e~^e~ J/ij) + gg at 
0{al). It turns out that this correction is comparable in magnitude with the NLO QCD 
correction, if not more important. 



^ The end-point coUinear logarithms in the color-singlet channel has been resummed, but its impact on the 

J/tp spectrum seems rather insignificant [45, 46]. 

However, although including NLO QCD correction helps to get the right answer for the inclusive produc- 
tion rate and energy spectrum of J/tp in e+e" — )■ J/ip + Xcc, it was noted that even [44], it is still difficult 
to reproduce the measured J/tp polarization and angular distribution. 



^ The relativistic correction to e+e J/tp + Xcc has been calculated and was reported to be surprisingly 
small [22]. 



a[e+e- ^ J/^ + X.g] = 0.74 ± 0.08^[;:[;^ pb, 
(T[e+e- ^ J/ip + Xiight] = 0.43 ± 0.09 ± 0.09 pb. 



(la) 
(lb) 
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An experienced reader may agree that, calculations of QCD perturbative corrections 
can be guided by some standard and unambiguous procedure. By contrast, calculating 
relativistic corrections, unexaggeratedly speaking, seems often plagued with ambiguities 
and pitfalls ^. The problem gets particularly acute, when the kinematics becomes involved, 
as in our case with three-body final states. It often occurs that, different results have been 
reported by different authors in calculating relativistic correction to the same process ^. 
To the best of our knowledge, a simple and consistent recipe for calculating relativistic 
corrections for a generic quarkonium production process has not yet been explicitly given 
in hterature. One of the major motifs of this work is also to fill this gap. We attempt to 
utilize a convenient yet slightly unconventional strategy to deduce the relativistic corrections, 
applicable to any inclusive quarkonium production (decay) process in color-singlet channel. 
Our approach is slightly different from the orthodox NRQCD doctrine, in that the matching 
coefficients are not truly of "short-distance" nature, for they explicitly depend upon physical 
kinematic variables rather than the partonic ones. However, we stress our method is still 
consistent, for its validity resting upon a rigorous relation in NRQCD, the Gremm-Kapustin 
relation [53]. 

The rest of the paper is structured as follows. In section II, we state the NRQCD ex- 
pansion formula relevant to this work and discuss the corresponding long-distance matrix 

elements of the color-singlet production operators. In section III, we outline our matching 
scheme that can be applied to any color-singlet quarkonium production process, and dis- 
cuss its advantage over the more orthodox doctrine. In section IV, we give the differential 
expressions for the three-body phase space needed for the reaction e'^e~ J/ip + gg. In 
section V, we present a detailed description on how to determine the short- distance coeffi- 
cients through relative order v"^ in e+e^ J/i'' + Xnght and how the physical predictions 
for the inclusive production rate of J /ip in e^e^ annihilation come out. In section VI, we 
apply our formulas to investigate the phenomenological impact of the order-t»^ correction to 
the integrated production rate and the differential energy spectrum for the unpolarized J/V', 
and the energy distribution of angular and polarization parameters at B factories. Finally 
we summarize our results in section VII. In Appendix A, we collect the analytic expressions 
for numerous types of differential cross sections of J /ip in e^e~ annihilation, at the leading 
order and next-to-leading order in v"^. In Appendix B, we show that, it is possible to reex- 
press our predictions for the integrated cross sections in a more orthodox form, that 
depending explicitly on the charm quark mass rather than the J/ if) mass. 

II. NRQCD FACTORIZATION AND LONG-DISTANCE MATRIX ELEMENTS 

According to the NRQCD factorization, the inclusive J/^ production rate in e+e" colli- 
sion can be schematically written as 

(ia[e+e- ^ J/V^ + AT] - ^ da[e+e- cc{n) + X]{Oi''''). (2) 

n 



^ See Ref. [47] for an early discussion on one type of ambiguity affiliated with the normalization of particle 

states in relativistic correction calculations. See also [48] for a related discussion. 
^ For example, Ref. [49] claims a disagreement with an earlier publication [50] on the result of relativistic 
correction to photoproduction of J /ip- Likewise, a recent calculation of the relativistic corrections to the 
fragmentation function for the c quark to fragment into J/V' [51] also disagrees with an earlier work [52]. 
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The NRQCD expansion is organized by the velocity scahng of the vacuum matrix element 
of NRQCD 4-fermion operators, {OiJ'^)^ where n denotes the collective quantum numbers of 
the cc pair created in the hard scattering. The ddn are canonically referred to as the process- 
dependent short-distance coefficients, which depend on all partonic kinematic variables for a 
given production process, but are insensitive to the long-distance aspects of the quarkonium 
state J/ip. 

In this work we are only concerned with the color-singlet channel. The readers who are 
also interested in the color-octet contributions to this process can refer to Ref. [36, 37, 39]. 
For our purpose, the relevant 4-fermion color-singlet production operators are given by ^ 



1 



2ml (2J+ 1) 



(3b) 



where -0 and x are Pauli spinor fields for annihilating and a heav^ quark, and creating a 
heavy antiquark, respectively, a* denotes the Pauli matrix, and D is the spatial part of 
the antisymmetrical covariant derivative: V'^Dx = '^/'^Dx — (DipYx^ ^ form to preserve 
Galilean invariance. J — 1 denotes the total spin of J/ijj, and X denotes additional fight 
hadrons accompanied with J/ip with net energy no larger than the ultraviolet cutoff of 
NRQCD. Note the sum for the intermediate states is extended not only over all additional 
light flavor states X, but over 2 J + 1 polarizations of the J/ip as well. 

The vacuum expectation values of these NRQCD production operators are genuinely 
nonperturbative objects, whose exact values are even difficult to ascertain from the powerful 
nonperturbativc tools such as lattice QCD, mainly owing to the obstacle in implementing 
those asymptotic states containing X. Fortunately, in practice one can always invoke the 
so-called vacuum saturation approximation (VSA) for the color-singlet channel, which is 
accurate up to an error of relative order to link these NRQCD operator matrix elements 
with the more familiar Schrodinger wave functions at the origin for the J/ip: 

{Oi'^) ^ (J/^|Oi(^5i)bbl|JM ^ |(0|xWl^/^)r = 2iVe^}/^(0), (4a) 



{Vi'"^) « (J/V'|Pi(^5i)bbl| JM « Re (J/V^lV^VxlO) • (0|xV(-|D)>| J/V') (4b) 



-2N^ Re 



V'}/^(0)VVj/^(0) 



Under this approximation, the vacuum matrix clement for J/ if) production can be approx- 
imated by the square of vacuum-to- J/t/^ matrix element in NRQCD, and further by the 
corresponding decay matrix element for the J/%1) state. 

The leading-order color-singlet J/'^ production matrix element, (Oi^^), can be identi- 
fied with the familiar wave function at the origin for J/ip, ipj/^{0). This quantity can be 



^ In this work, wc find it convenient to choose a different normaiization for the J/xp production operators 
0'(^^ and P'^^^' other than the standard ones introduced by Bodwin, Braaten and Lepage (BBL) [3], i.e., 
we define o(^'' = 2jTT*-^/bbl ^^'^ 'P'^/'^ = (2 j+i) ^/bbl ■ Note the prefactor 1/m^ normalizes the 
operator p^^'^ such as to carry the same mass dimension as O^^^ . 
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determined by several means, e.g., from lattice simulation, or from plienomenological quark 
potential models, or directly from the measured leptonic width of J 

The determination of the relative order-f^ production matrix element, (P/^^), turns to 
be more problematic. In Coulomb gauge, the gauge field piece in this matrix element is 
suppressed relative to the ordinary derivative piece. By VSA, it seems intuitive to interpret 
{Vi^^) as product of V'J/i/'(0) ^^"^ the second derivative of the wave function at the origin, 
V^'?/'j/^(0). Nevertheless, one should be cautioned that such a naive interpretation is obscure. 
This is because the bare NRQCD matrix element contains linear ultraviolet divergence, hence 
it needs to be regularized and renormahzed, thus depending on the cutoff of the NRQCD 
lagrangian (An overbar put above the wave function is to remind this). There seems no 
direct way to directly infer this matrix element from phenomenological potential model. 
Nevertheless, it is the NRQCD effective theory framework that endows this nonperturbative 
quantity a meaningful definition. 

For later use, it is convenient to introduce the dimensionless ratio of the vacuum matrix 
elements of the following NRQCD operators: 

, 2. ^ (J/^(A)|^t(_l6)2^.,(A)^|0) 

^^^^/^ ml{J/^l^iX)\^^a■eiX)x\0) ' 

This quantity, characterizing the typical size of relativistic correction for J/ip, is supposedly 
around 0.3. Note that its value is independent of the J/ip helicity A in above equation. 

The vacuum-to-quarkonium relativistic correction matrix element has been measured by 
lattice QCD, though the uncertainty is quite large. There is an interesting relation, first 
derived by Gremm and Kapustin (G-K) [53], which derives from the equation of motion of 
NRQCD, and expresses the relativistic correction NRQCD matrix element in terms of the 
LO NRQCD matrix element, physical J/ip mass and the charm quark mass: 

^ - l + l{v')j/, + 0{v% (6) 

In NRQCD, the quark mass parameter is most naturally identified with the quark pole mass. 
Unfortunately, due to the intrinsic ambiguity of the charm quark pole mass, this relation 
cannot be utilized to nail down the precise value of ('y^)j/^. It can not be precluded that, 
the actual value of this quantity might be quite far from the naive expectation, 0.3. Since it 
is a subtracted quantity, it will be perfectly consistent if {v^)j/^ turns to vanish or become 
even negative ^. 

During the era anteceding the inception of the NRQCD approach, many authors preferred 
to using the binding energy, i.e. e = Mj/,p — 2mc to parameterize the contribution of 
relativistic corrections (for example, see [50, 55]). With the aid of the G-K relation (6), all 
those old results can be readily translated into the modern form, with relativistic correction 
designated by the NRQCD operator matrix elements. 



^ It is worth mentioning that, recently there have been claims that this quantity can be quite accurately 
extracted from the Cornell-potential-model-based analysis, = 0.225lQ;Qgg [54]. 
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III. PERTURBATIVE MATCHING STRATEGY 



The central ingredient of the NRQCD factorization formula is to deduce the NRQCD 
short-distance coefficients. The procedure of determining these coefficients are usually re- 
ferred to as matching. The idea is rather straightforward, since these short-distance coef- 
ficients are in principle insensitive to the long-distance nonperturbative physics, therefore 
one may replace J/ip hj a, free cc pair carrying the quantum number of '^S^\ then both 
sides of Eq. (2), including the NRQCD matrix elements in the right side, can be accessed 
entirely in perturbation theory. Matching both sides, one then be able to extract the desired 
short-distance coefficients dcFn- 

Let us take e^e~ J / ip + X a,s an explicit example to illustrate the problem faced for 
the matching calculation beyond the lowest order in v"^. Schematically, one can express the 
corresponding perturbative matching formula for (2) as 

Y,{'^T^Y5\K-P-kx)\M [e+e-^ccf>5f\P,A)+x]|' = J] da„(P, A)(e'5, (7) 

X n 

where the flux factor associated with the single-inclusive cross section has been suppressed 
for simplicity. K stands for the sum of momenta of the colliding electron and positron, i.e., 
the 4-momentum of the virtual photon into which the electron and the positron annihilate. 
The sum in the left side is extended over the spins of all the additional partonic states X, 
as well as over the phase space integration affiliated with X. 

It is clear from (7) that, to identify the matching coefficients, one needs to expand the 
inclusive production rate for perturbative cci^S'^'') pair systematically in the small relative 
momentum between c and c, q. This procedure entails two essential ingredients, one is 
to expand the matrix element squared in powers of q, the other is to expand the phase 
space integrals accordingly. The former operation is more or less standard, but the latter 
potentially cause some problems for the orthodox matching method, as will be reviewed in 
section 111 A. The main trouble is that, in the standard matching calculation, it is mc, instead 
of the physical J/%Ij mass, Mjj^ that should enter the NRQCD short-distance coefficients. 
The orthodox method then requires that the phase space integral be expanded around a 
fictitious cci^S^'^) state of mass 2mc. Technically, such an expansion of the phase space 
integral at differential level is not easy to realize. More importantly, this operation leads 
to some inevitably unsatisfactory feature in predicting the differential J/ip spectrum: when 
approximating J/%Ij mass by 2mc, the incorrect kinematics causes the spectrum somewhat 
distorted, which may become particularly problematic near the phase space boundary 

In section IIIB, we will elaborate on the matching method adopted in this work. Mo- 
tivated by the aforementioned shortcoming of the orthodox matching strategy, we attempt 
to circumvent the most difficult part arising from expanding the phase space integral. The 
key point is that we choose to use physical kinematics instead of the partonic one, and the 



If one is content to knowing only the total cross section, this orthodox method should be straightforward 
and does not cause any problem. For instance, for simpler reactions such as gg ^ r]c, exclusive double 

charmonium production e+e~ — > J/ijj + rjc, or inclusive quarkonium decays J/ijj — > ggg, it is possible to 
first work out the phase space integration analytically, then expanding the resulting integrated partonic 
production/decay rates in powers of q about a fictitious charmonium of mass 2mc. However, in the case 
of more involved kinematics, it is usually not feasible to acquire the integrated rate in a closed form. 
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invariant mass of the cc pair appearing in the matching calculation is taken as Mj/^. As 
we shall see, this brings in great technical simplifications. As a result, we can perform the 
matching at the level of the matrix element squared, instead of at the level of the produc- 
tion rate. Although our matching method somewhat deviates from the ordinary tenet in 
that the "short- distance" coefficients now explicitly depend on Mj/^, it is nevertheless still 
theoretically consistent, thanks to the G-K relation (6). 



A. Orthodox matching strategy motivating the shape-function method 

In this subsection we review what the standard NRQCD matching strategy would look 
like. Historically, this method antecedes, and, motivates, the so-called shape function 
method [15]. The orthodox doctrine of NRQCD matching is common in any effective field 
theory, in that the short-distance coefficients should depend only on the parton kinematics, 
thus on the quark mass, and there is no way for quarkonium mass, which inevitably entails 
the long-distance hadronization effect, to enter into them. 

Let c and c that evolve to J/ip in (7) have momenta p and p. Both c and c are supposed 
to be on-shell, and their momenta can be decomposed as 

p = ^P + qi, (8a) 
P = Ip + Q2. (8b) 

Here the "total" momentum P, which is dehberately chosen to satisfy a/p^ — 2mc, should 
be distinguished from the true total momentum of the pair, P = p -\- p, with invariant mass 
of 2Eq, where Eq = \J mf. + is the energy of the c or the c in the rest frame of the cc 

pair. In the rest frame, these 4- momenta have following explicit assignments: P'^ = (2mc, 0), 
Qi = {Eq — rric, q), q'2 = {Eg — mc, — q), respectively. In the laboratory frame, it is understood 
that a suitable boost along the moving direction of the cc pair is imposed on all the 4-vectors. 

The purpose of introducing P is that one needs to expand the phase space integral around 
a fictitious cc pair of invariant mass 2mc, which serves as the basis momentum. Concretely, 
the constrained phase space measure for the partonic process of (7) is 

where K is the momentum of the virtual photon, and ki represents the additional partons 
in X. Note Qi inside the 5-function are understood to be subject to an appropriate Lorentz 
boost from the rest frame of P to the laboratory frame. 

The squared quark amplitude can then be folded with the phase space measure (9) to 
obtain the partonic cross section. The cross section needs to be expanded in the small 
momenta q^, and powers of momentum can be identified with derivatives acting on the 
heavy quark fields according to NRQCD factorization. Factors of relative momentum qi — q2 
(in the rest frame of the cc pair) typically arise from expanding the amplitude, which can 

be identified with the 'ip'{i'D)x- Furthermore, in expanding the 5-function in phase space 
measure (9), one typically encounters a different type of factor, the center-of-mass (cms) 
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type momentum g° + Q'2 i^^ ^^e rest frame of the cc pair), which can be identified with a 
total time derivative iDo('^^x)- 

As noted in Ref. [15] (see also [16]), these cms-derivative operators, though nominally of 
higher-order than the relative momentum operators in NRQCD expansion, can be of dynam- 
ical significance near the kinematic endpoint of quarkonium spectrum. Upon expansion of 
the 5-function in qi + q2, the resulting power series in make increasingly singular contribu- 
tions near the boundary of partonic phase space, which signals that NRQCD expansion may 
break down near the endpoint region. Fortunately, it has been shown that such enhanced 
kinematic contribution due to these cms relativistic corrections can be resummed, whose 
effects are then encoded in the universal nonperturbative shape function [15]. 

The shape function is of greatest utility to improve the predictions for inclusive quarko- 
nium production in the color-octet channel [15, 16, 39]. Nevertheless, it can also play a 
nontrivial role even for the color-singlet channel, which is relevant to our case. It turns 
out that the resulting series from expanding the 5-function in (9) can be exactly resummed 
without introducing any new nonperturbative parameter other than the quarkonium mass. 
Its sole effect is to account for the difference between quark and quarkonium mass, conse- 
quently shift the partonic boundary of phase space to the hadronic one. The remarkable 
effect can be easily understood. The cms-momentum factor + equals 2Eq — 2mc in the 
rest frame of P. When identified with the total time derivative i9o('?/'^x), this operator can 
convert to the binding energy e = Mj/^ — 2mc when sandwiched between the vacuum and 
the physical J/ip states, thus helping to recover the hadronic kinematics. Not surprisingly, 
the underlying reason is nothing but the G-K relation. 

The role played by the color-singlet shape function seems to strongly suggest that, the 
inconvenience brought in by the orthodox matching method, i.e.. the procedure of expanding 
and reassembling of the phase-space 5-function, may be avoidable. As will be elaborated 
in more detail in next subsection, one may just remain the physical kinematics intact in 
(9) throughout the matching computation. Lastly, we mention the fact that, somewhat 
ironically, there seems no any practical calculation of the complete first-order relativistic 
correction for the inclusive quarkonium production process that is based on this orthodox 
matching tenet. 

B. Matching strategy adopted in this work 

In light of the complication inherent in the orthodox matching method, in this section 
we are going to describe a different matching strategy, which is suitable for any inclusive 
quarkonium production (decay) process in the color-singlet channel. We will take the reac- 
tion e~ e'^ ^ J/ip + gg as an explicit example to illustrate our method. 

When assigning the momenta of c and c in perturbative matching, the separation between 
"total" and "relative" momenta is just a matter of taste, by no means unique. Here we will 
choose a different one from that in section III A. The momenta of the on-shell c and c can 
be decomposed in the following form: 



We stress here P represents the true total momentum of the pair, P — p + p, with invariant 



P 



P 




(10b) 



(10a) 
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mass of 2Eq, and now P and q are chosen to be orthogonal: P • g = 0, in contrast to 
the choice made in (8). In the rest frame of the cc pair, the exphcit components of the 
momenta are = {2Eq, 0), = (0, q), = {Eg, q), and = {Eg, — q), respectively. In 
the laboratory frame, one has P^ — (P°,P) = (y^P^~+4E|, P) and appropriate Lorentz 
boost is understood to be imposed on any other 4- vector. It is worth mentioning that, it is 
this form of momentum assignment, rather than (8), that has been practically used in most 
calculations of relativistic corrections [22, 51, 56-59]. 

The greatest advantage of this kind of momentum decomposition is that, there is no need 
to expand the total momentum P of the cc pair around a basis momentum with invariant 
mass of 2mc, and we will leave phase space measure intact by assuming the cc pair with an 
invariant mass 2Eg. We argue by this way the relativistic effects in phase space measure are 
automatically incorporated, at least to relative order v^. We will come back to the connection 
between the factor Eg and Mj/,^ in nonrelativistic expansion later in this subsection. 

Even though we are coping with inclusive J/ip production process, but insofar as the 
color-singlct channel is concerned, it is not necessarily be committed to the cross section 
level at the very beginning. In fact, since we no longer need worry about the complication 
from the phase space integral, it seems legitimate to invoke the NRQCD factorization at 
the amplitude level. To this end, we need only retain those operator matrix elements that 
connect the vacuum to the color-singlet J/ip state. To the order of desired accuracy, the 
amplitude can be written as 



M[J/iP{X) + gg] = y^2Mjj;[Fo{X){J/^\^''a-ex\0) 

(11) 



+ ^WIV'V-6(-|Dfx|0) + - 



where -Fi(A) are the corresponding short-distance coefficients, which are Lorentz scalars 
formed by various kinematic invariants in the reaction. In particular, they also depend 
explicitly on the helicity of J/V", A. For the Lorentz- invariant amplitude in the left-hand 
side of Eq. (11), Ai[J/ip{X) + gg], it is most natural to assume relativistic normalization 
for each particle state, since the squared amplitude needs to be folded with the relativistic 
phase space integral to obtain the physical cross section. However, in the right-hand side of 
Eq. (11), the J/ip state appearing in the NRQCD matrix elements conventionally assumes 
the nonrelativistic normalization. To compensate this difference, one must insert a factor 
y'2Mj/^ in the right side of (11). 

Squaring both sides of (11), summing over the final-state spin/colors as well as averaging 
upon the initial-state spins, the matrix element squared reads 

Y,\M[J/^P{X) + gg]\' = 2Mj/40(/'^)Y,{\Fo\' + 2Re[FoF;]{v')j/^ + ---}, (12) 

where the VSA has been invoked and (Of^^) has been given in (4a), and {v'^)j/^ defined in 
(5). The symbol ^ indicates the suitable color-spin summation/average. 

To determine the coefficients |FoP and Fo-^2* + H-C-, we follow the moral that these short- 
distance coefficients are insensitive to the long-distance confinement effect, so one can replace 
the physical J/ip state by a free cc pair of quantum number ^S[^\ by which the NRQCD 
operator matrix elements can be perturbatively calculated. The short-distance coefficients 
Fi{X) can then be read off by comparing the full QCD amplitude for producing cc{^S^^) and 
the corresponding NRQCD factorization formula. In our case, the amplitude for producing 
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a cc(^S[^'') pair associated with two gluons is 

M[ccCsi'\ A) + gg] = FoiX){ccCS^)\i;^a ■ ex\0) + ^{ccCSM^a ■ e{-^3fx\0) 



mi: 



Fo(A) + F2(A) 



(13) 



In Eq. (13), we use relativistic normalization for the c and c states in the computation of the 
full QCD amplitude and in the computation of the NRQCD matrix elements. Consequently, 
a factor 2Eq appears in the second expression of Eq. (13). An additional factor \/2Nc 
arises from the spin and color factors of the NRQCD matrix elements. Prom (13), it is 
straightforward to extract the short-distance coefficients Fi[\): 



F2iX) 



M[cc{'Si'\X) + gg] 



ml (M[cx:('S^^\\)+gg] 



V2K2Eg 



-Fo{X) 



(14a) 
(14b) 



q2=0 



The LO coefficient Fq can be obtained by putting q — )■ in the amplitude and equating 
Eg and rric. To deduce the coefficient F2, it is understood that one has to first expand 
the amplitude to the first order in prior to taking the q ^ limit. Consequently, it is 
necessary to distinguish between Eg and rric- 

Although the expression of Fq can be unequivocally determined, it is not without ambigu- 
ity to deduce the coefficient F2. This is because, determination of this relativistic correction 
coefficient crucially hinges on which quantity is chosen to be expanded around in powers of 
q^ in the quark amplitude. 

Obviously, those terms that contain exphcit q^ factor should contribute to F2. Besides 
these terms, in the matching procedure adopted by most authors, one usually often includes 
relativistic effects implicit in all the expressions that contain the factor Eg, where Eg is 
always expanded around rric in power series of q2, i.e. Eg = mc + ^^+ 0{q*). By collecting 
all the sources proportional to q^, one is then able to deduce the coefficient F2 according to 
(14b). 

In this work, we find it much more advantageous to take a somewhat different route. 
Aside from retaining those relativistic correction contributions that contain q^ explicitly, we 
choose to expand every occurrence of rric in the amplitude in term of q^/E^, while keeping 



Eg intact: 



m, 



C — Eg 



2Er, 



+ 0(q')- 



(15) 



Somewhat nonstandard as it seems, but as we will see shortly, by choosing this way of 
expansion, we circumvent the most difficult task, i.e., expanding the three-body phase space 
integral. This procedure turns out to be the simplest in practice, especially when contrasted 



Throughout this work, the bold-faced symbols, such as q, if not otherwise stated, arc exclusively referring 
to the spatial vectors defined in the rest frame of the cc{P) pair, whereas the italic symbols, such as q, 
are reserved for the covariant 4-vectors, often presumed in the laboratory frame. 
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with the orthodox matching method outhned in section III A. This will be exemplified by 
more concrete examples in section V. 

In (14b), we have refrained from expressing F2 by taking the second-order derivatives 
of the quark amplitude over as frequently adopted in many works [22, 51, 56, 58]. The 
reason is that we try to avoid potential ambiguity associated with this operation, which 
usually happens when one performs the standard expansion around itLc. The recipe given in 
(14b) is unambiguous and simple provided that Eq is kept fixed. The expression obtained 
this way are connected to the standard one through reshuffling some terms between Fq and 



Squaring the matrix element (13), summing over final-state polarizations and averaging 
upon the initial-state spins, we obtain 



The matrix elements (Cf^) and {V^'^) denote the vacuum expectation values of the produc- 



wherc the factor of 2Nf. in the right side of Eq. (17a) arises from the spin and color normal- 
ization factors for free cc states. Comparing both sides of (16), one may directly deduce the 
short- distance coefficients \Fq\^ and Re[FoF|]. 

In passing it may be worth reminding that, during the polarization sum/average proce- 
dure, new factors of Eq will be unavoidably regenerated in the squared amplitude. Evidently, 
such factors can arise from summing the polarization states of cci^S^^) state. In the stan- 
dard way of expansion, one needs re-expand these occurring Eq factors once and more, and 
keeping reshuffling the corresponding terms from the LO matrix element squared to the 
relativistic correction piece. Fortunately, since we keep Eq fixed in our approach, no any 
extra labor needs to be invested for such complication. This comprises another attractive 
trait of our expansion strategy. 

Substituting the short-distance coefficients (14) to (12), or directly converting the quark 
amplitude squared (16) to (12), after some straightforward algebra, one then obtains the 
desired squared matrix element for producing J/%1) plus light hadrons. 

There arises one immediate question. Since rric has been eliminated in favor of Eq in 
the physical matrix element squared, it is necessary to specify which value of Eq should be 
taken, in order to make concrete predictions. If there were no rationale to restrict its value, 
our approach would just yield ad hoc predictions and lack attractiveness. 

Fortunately, the answer to this question is definite, i.e., theoretical consistency requires 
that Eq can be fixed in an unambiguous manner. To see this, let us first make the following 
observation: 



F2. 



J2 \M[cc{'S„X)+gg]\' = AE'^Yl {\m)\'{Of) + 2Re[FoF;]{Vf) + ■ ■ ■} 

= 4£;^2(27V,)^||Fo(A)r + 2Re[FoF,*]^ + ---|. (16) 




(17b) 



(17a) 




(18b) 



(18a) 
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where (18a) comes from simple nomelativistic kinematics, the inverse relation of (15). 
Eq. (18b) is nothing but the G-K relation (6). 

The similarity between (18a) and (18b) strongly suggests that, Eg appearing everywhere 
in the short-distance coefficients in (12) can be replaced by Mj/^/2. Naively, the entering 
of J/ih mass into short-distance coefficients seems to be a nuisance, which is against the 
doctrine of the EFT. Nevertheless, this procedure is valid, at least to the present accuracy 
of order v"^, thanks to the G-K relation 

Identification of 2Eq with Mj/^, in conjunction with our nonstandard expansion (15), turn 
out to have great advantages. By this way, the relativistic effects in phase space integrals 
are automatically incorporated. In some sense, our approach fulfills the role of the color- 
singlet shape-function by promoting the partonic kinematics to hadronic kinematics, but 
not necessarily restricted to the region of the maximum J/tjj energy. In addition, since the 
mass of J/ip is known rather precisely, it is better to choose it as the input parameter than 
the ambiguously defined charm quark mass. To summarize, our method greatly simplifies 
the efforts required for the matching calculation, by reducing the task of matching the cross 
section to matching the amplitude squared. 



IV. THREE-BODY PHASE SPACE FOR e+e" cc{^Si) + gg 

One integral part of the matching procedure is to consistently incorporate the relativistic 
effects inherent in the phase space integration. As was explained in section IIIB, owing to 
the virtue of our matching approach, no special care needs to be paid to the phase space 
integral, provided that we identify the invariant mass of cc pair, 2Eq, to be Mjj^. For the 

process e~(/i) -|- e+(/2) — 1*{K) — )■ cci^Sf^ ^ P) -|- g{ki) + g{k2) considered in this work, the 
energy- momentum conservation requires K — li -\- 1'^ — P -\- ki -\- k^- The electron and gluon 
are treated to be massless, and = ^'j/ti>- ^^^^ evaluate the three-body phase space 
(ilia in the e+e" center-of-momentum (laboratory) frame. 

The center-of-mass energy squared is defined by s = . It is also convenient to define a 
dimensionless ratio r = M'^i^/s. The differential three-body phase space can be expressed 
as follows: 

l+r pi rxf 1-2-K 

dz d cose I dxi / d^l. (19) 

2^ J-1 Jx^ Jo 



2(47r) 

It is worth noting that, in contrast to (9), the advantage of our matching method is there is 
no need to expand P around the momentum of a fictitious particle with mass of 2mc- For 
notational simplicity, we have introduced three dimensionless variables, z = = 

Xi = ^2 = X2 = ^2 = to characterize the fractional energies carried by the 

cc{^Sl^^) pair, the gluon 1, and gluon 2 in the laboratory frame. Only two of these three 
variables are independent, since they are subject to the constraint from energy conservation: 



It is worth noting that this substitution has also been adopted in a recent investigation of the relativistic 
corrections to the exclusive charmonium production process e+e~ J/tjj + r]c [23]. 
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xi + X2 + z = 2. In the following, we will eliminate X2 everywhere in favor of z and xi as 
the integration variables. 

We use {9, 4>) to denote the polar and azimuthal angles of the outgoing J/ip momentum 
with respect to the moving direction of e~ in the laboratory frame. We have suppressed d(f) in 
the integration measure since it has been trivially integrated over due to the axial symmetry 
of the reaction under consideration. It is convenient to introduce a set of auxiliary solid-angle 
variables (0^, as the polar and azimuthal angles of the moving direction of the gluon 
1 in a rotated coordinate system relative to the laboratory frame, where the J/ip moves 
along the new +z axis. Given the energy fractions z and xi, four- momentum conservation 
uniquely constrains the polar angle O^: 

2(l + r-.)-..(2-.) 

XiV z'^ — 

The main advantage of choosing these integration variables as given in (19), is that each 
of them has an intuitive interpretation and the respective integration boundaries are rather 
simple. This is in contrast with the set of variables employed in Ref. [22, 28, 59]. 

The integration boundaries for z have been explicitly labeled in (19), and those for xi 
can be easily inferred: 

. ^-^f^ . (21) 

For the suppressed variable a;2, the boundaries would be the exactly same as x^. 

If one concentrates only on the energy spectrum of (un)polarized J /'ip, disregarding its 
angular distribution, one may take a shortcut- by starting with the simpler process 7* — )■ 
J/'0 -|- gg, then converting the differential decay rate to the J /ip differential production 
cross section [25]. In such case, since there is no preferred orientation in space, two angular 
variables, cos^ and $1, can be trivially integrated over in (19), consequently one is left with 
only two dimensionless energy variables in the three-body-phase-space measure: 



2(47r) 



l+r fX 



dHa = / dz / dx^. (22) 



In some situation, it is desirable to know the analytic expression for the integrated J/-?/' 
production rate. To this purpose, it seems more advantageous to choose a different order to 
perform the phase-space integration: 



l+r 

rfHa = ;^7^ / dxx I dz. (23) 



2(47r)= 



As an intermediate byproduct, one can deduce the gluon energy spectrum once performing 
the integration over z. Phenomenologically, knowing this is not so meaningful, since it 
cannot be directly linked with a physical observable. However, as a calculational device, 
choosing this particular order for phase-space integration leads to considerable technical 
simplicity, because the integration boundaries in (23) are far simpler than those in (22). As 
a consequence, by this way one can readily deduce the J/ ip total production rate in a closed 
form, which is otherwise rather difficult to achieve if one starts from (22). 
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V. OUTLINE OF MATCHING CALCULATIONS FOR e+e" cc{^Si) + gg 



In this section, we present a detailed description on how the short- distance coefficients 
through relative order-f ^ can be determined via our matching procedure, concretizing the 
method put forward in section IIIB. We also illustrate how the various types of predictions 
for the inclusive J/ip production rate in e~^e~ annihilation emerge. 

In order to deduce the intended short- distance coefficients, one needs to consider the 
parton process e+e~ — )■ cc{^S[^\P,X) + gg, with one typical lowest-order diagram shown 
in Fig. 1. The calculation is expedited by the covariant projection technique developed by 
Bodwin and Petrelli [57], which helps to readily project out the amplitude for the cc pair 
being in the color-singlet spin-triplet state 

Our matching method will be exemplified by the following three subsections. In sec- 
tion VA and section VB, where we are only interested in the energy spectra of unpolarized 
and longitudinally-polarized J/ip, we take the shortcut by considering the simpler process 
7* — )■ cc{^S[^\ P, A) + gg; in section VC, where we are also interested in the angular as well 
as the energy distributions of J/ip, we work with the full process e~^e~ — )■ cc{^S[^\ P, X)+gg. 




FIG. 1: Lowest-order Feynman diagrams for e^e — >• J/ip + gg- 



A. Matching calculation for 7* cc{^Si) -\- gg 

We start with the tree-level quark amplitude J*{K) — > c{p)c{p) + g{ki) + (7(^2)) with the 
momenta of c and c defined in (10) The amplitude can be written as 

u{p)s^v{p) = Tr [v{p)u{p)j^] . (24) 



An alternative approach, the threshold expansion method [60], is also valid to deduce the NRQCD matching 
coefficients. This method has been utilized to investigate the order-w^ relativistic correction to J/ip 
photoproduction at HERA [49]. However, for the process with involved kinematics like ours, this method, 
which requires to extensively deal with the algebra of two-component spinors, seems not as convenient as 
the covariant projection technique outlined in [57]. 

The calculation presented in this subsection is somewhat similar to that for the process g* — >■ J/ip + gg, 
from which the relativistic correction to the gluon-to- J/?/' color-singlet fragmentation function can be 
extracted [56]. 
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Here is a Dirac-color-space matrix, which reads 



^ = {e.eg'^ T'^T' ® f{g; h)-/-— ^(7*; K)^- f{g- h) + 5 perms, (25) 

|) + f 1 - mc -p- f2-mc 

where Qg signifies the QCD couphng strength, and eCc denotes the electric charge of the 
charm quark (cc = §)■ a, 6 denote the color indices of the two gluons, 6(7*; X), e*{g;ki) 
and e*{g;k2) represent the polarization vectors of the decaying virtual photon, gluon 1 and 
2, respectively. 

One of the important sources of relativistic corrections stem from expanding the quark 
propagators. Apart from retaining the factor q in the numerator of the propagator, one 
needs also expand its denominator to the quadratic order in q. Taking the diagram shown 
in Fig. 1 as an example, two propagators there are expanded to be 



1 2q-ki 4(g • ki) 



2 



+ T^^ + 0{q% (26a 



(p + kiy-ml P-ki {P-kiY {P-kiY 



(p + k2y-ml P-k2 {P-k2Y {P-k2f 

To proceed, we need project the amplitude (24) onto the spin-triplet color-singlet c{p)c{p) 
state, by replacing the v{p)u{p) with a suitable projection matrix. The projector that is 

valid to all orders in q for the spin-triplet color-singlet channel, denoted by A.^^\p,p, A) (A 
characterizes the polarization of this spin-triplet cc pair), assumes the particular form [57]: 

4\P.P. A) = -7^777^7- - ^c) r (A)(f +2E,)(|^ + me) ® (27) 

where Ic is the unit matrix in the fundamental representation of the color SU (3) group, and 
the spin-polarization vector e*(A) satisfies P ■ e*(A) = 0. The above spin projector is derived 
by assuming the relativistic normalization convention for Dirac spinor: u^'^^u^^^ — 2mSrs and 
^^('•)t^(s) — 2Eq5rs- Applying this spin-color projector to (24), we obtain 

McciP,q,X;ki,k2) = TT{s^Af\p,P,X)}, (28) 

where the trace acts on both Dirac and color spaces. Aicc{P, q, A; ki, ^2) can be interpreted 
as the amplitude for producing a color-singlet, spin-triplet cc pair in association with two 
gluons. 

We have emphasized that our method differs from the conventional tenet of matching. 
Rather than Taylor-expanding Eg around m^ in q^/m^ everywhere in Aicc{P, q, A; ki, k2), we 
choose to expand rric around Eg in powers of q^/E^ using (15). It is worth reminding that 
we should not forget to trade rric for Eg that appears in the denominator of the projector 
(27): 



1 1 / 1 q 



.2 



It is convenient, at this stage, to truncate the amplitude A4cc such that all the terms in it 
are at most quadratic in q. 

In the amplitude (28), the cc pair is warranted to be in the spin-triplet, but not necessarily 
in the S-wecve orbital-angular-momentum state. To project out the S-wave amplitude, one 
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needs average the amplitude Aicc over all the direction of the relative momentum q in the 
rest frame of cc{P) pair. This literal angular averaging procedure can help to acquire a 
specific class of relativistic corrections to all orders in v. However it is feasible only for 
few processes with very simple kinematics [23, 57]. For the process at hand, this procedure 
would become extremely cumbersome, if not impossible. Fortunately, to the intended 0{v'^) 
accuracy, one can utilize a standard trick to project out the S'-wave part. Now we already 
have the amplitude truncated up to two powers of q. Terms that contain no powers of q, or 
contain exphcitly the Lorentz scalar (which can be translated into — q^), already yield a 
pure S'-wave contribution; for those terms containing the tensor q'^q'^ contracted with other 
4- vectors, we can make the following substitution to extract the S'-wave piece 

^ yn'^'^(P), (30) 

where 

W^{P) = -gt^^ + ^^. (31) 

Following all these steps, we then obtain the desired ^Si piece of the amplitude, 
A^(^S'i, P, A; fci, /C2), accurate through the order q^. Following what has been elaborated 
in section IIIB, at this stage it is legitimate to replace Eq everywhere by in 
M{^S^,P,\-hM)- 

It is now the time to deduce the desired "short-distance" coefficients Pq and P2, following 
the recipes given in (14a) and (14b). After this is done, we need square these coefficients 
and perform the corresponding spin-color sum/average: 

2^Re[PoP2*] = ln,,,{K)U,,,{P)i-g^^,)i-g^^,)2Re[j^i;-'''''^:F;'^^^ (32b) 
where the amputated "short-distance" coefficients (i = 0, 2) are defined through 

Fi = J't '"^^,{1*: K) e;CS,; P, A) el{g; h) e^g; k,). (33) 

For simplicity, we have suppressed the color indices, implicitly contained in J^j, which reads 
ItiT^'T^) / yfWl. In Eqs. (32), we have summed over polarization and color of the cci^S^^) 
pair and two gluons, and averaged upon three spin states of the virtual photon (note the 
prcfactor |). Polarization sum for two gluon states has been taken into account by the metric 
tensors —Qaa' and —gisp'-i and that for the virtual photon and J/%Ij by the polarization tensor 
W^''{K) and Wp\P). 



Note this S'-wave projection operation will regenerate factors of Eq in the denominator of the amplitude. 
Obviously it will not cause any trouble for our strategy of matching. In contrast, in the orthodox matching 
recipe, one is forced to reexpand those terms containing these newly-generated Eq factors, and reshuffling 
the corresponding terms from the "leading-order" piece to the relativistic correction piece. This further 
exhibits the merit of our method. 
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According to equation (12), we then obtain the squared matrix element for 7* — )■ J/ip+gg. 
Including the 3-body phase space measure (22), we can obtain the differential decay rate for 
unpolarized J/ip: 

where we have included a statistical factor of ^ to account for the indistinguishability of 
two gluons. 

To convert the differential decay rate in Eq. (34) into J/ip production cross section, one 
can use the formula [25] 

da[e+e- J/ip + gg] ^ 4na dr[j* J/ip + gg] 
dzdxi dzdxi 

It is not difficult to integrate over the fractional energy of gluon 1, xi, with the integration 
boundaries specified in (21) to acquire the energy spectrum for unpolarized J/ip. 



B. Matching calculation for 7* — >^ cc{^Si,X — 0) + gg 



It is also interesting to know how the polarization information of J/ip varies with its 
energy. To accomplish this, in addition to the differential energy spectrum for unpolarized 
J/ip as given in section VA, we also need know that for the longitudinally-polarized J/ip. 

In this section we outline how the corresponding matching calculation is carried out. 

The "short-distance" coefficients Fq and F2 can be obtained following the the same pro- 
cedure as outlined in section VA. Nevertheless, the longitudinal polarization vector of J/xp 
can be explicitly substituted by 

which satisfies ■ P = and ■ = — 1- 

We then square these coefficients and perform the corresponding spin-color sum/average: 

E l^ol' = l^,AK)elfS,)e,,,CSi) {- 9 aa'){- 9 ') J^o'' J'o , (37a) 

2^Re [FoF*] = ^ n,,,(X) el^i'S,) e^A'Si) {-9aa'){-9pp') 2 Re[J-^^ . 

(37b) 

Here the amputated "short-distance" coefficients J-^ are the same as what appear in Eqs. (32). 
Obviously the sum over polarizations needs not act on the cc(^5'i, A = 0) state. 

Apart from replacing Eq everywhere with Mj/^/2, we also need substitute — and 

|P| = ^V-z^ — 4r for in above expressions. We then follow (12) to obtain the squared 

matrix element for 7* — )■ J/V'l + gg- With this expression at hand, one can use (34) to 
infer the differential decay rate from a virtual photon, subsequently use (35) to deduce the 
corresponding differential cross section for producing the longitudinal-polarized J/ip in e+e" 
annihilation. 
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C. Matching calculation for e+e cc{^Si) + gg 



In section VA and VB, we have resorted to a shortcut by considering the production 
rate oi a. J/ip + gg from a virtual photon decay, since the sole purpose is to deduce the 
energy spectrum of unpolarized or longitudinally-polarized J/ip in e~^e~ annihilation. In this 
subsection, we are interested in knowing the angular-energy double differential distribution of 
unpolarized J/ip. To this end, it is compulsory to begin with the full process e~(/i)e"'"(/2) 
J/tp{P) + g{ki)g{k2). We use h and I2 to signify the momenta of e~ and e"*", respectively, 
and I1 + I2 — K. 

The main result derived in section VA can be directly transplanted here. In particular, 
the "short-distance" coefficients Fq and F2 for 7* — )■ J/ip+gg, determined there by employing 
(14a) and (14b), only need undergo some slight modifications to meet our purpose. That is, 
one needs replace the polarization vector of the virtual photon by a e'^e" bispinor and insert 
a photon propagator and a QED coupling. These slight changes are embodied in squaring 
these coefficients and performing the corresponding spin-color sum/average: 



2^Re[FoF*] = 1 L,,,^ n,,,(P) (-w)(-^^/30 2 RefJ"^^ ''"^J^^ '''''^'''1 ■ (38b) 

Here the amputated "short-distance" coefficients J-^ are the same as what appear in 
Eqs. (32). The factor 1/4 represents the average over polarizations of the initial e~e^ 
state, and L'^'^' denotes the leptonic tensor: 

L^^^^' = E[^;(/2;r)7'^M(/i;s)][M(/i;s)7^'^;(/2;r)] 

s,r 

= A{l^,l^^ + 1^,1^; -h-hg^^'), (39) 



where the sum is extended over all possible polarization states of the electron and positron. 

Substituting Eqs. (38) into (12), we then obtain the color-spin averaged/summed matrix 
element squared for the process e+e^ J/ip+gg. Including the 3-body phase space measure 
(19) and the flux factor, we can obtain the differential production rate for unpolarized J/ip: 

In the squared matrix elements, all the scalar products can be expressed in terms of 
z, 6, Xi, and one additional angular variable, 61, which represents the angle between the 
3-momentum of gluon 1 and the beam direction in the laboratory frame. This polar angle 
is connected to 9, Ql and through 

cos 9i — cos 9 cos 01 — sin 9 sin cos , (41) 

where cos0^ is uniquely determined when z and Xi arc given, as indicated in (20). 

As first elaborated in [28], general consideration based on Lorentz invariance, parity and 
gauge invariance demands that for inclusive J/ijj production in e~^e~ annihilation, the double 
differential distribution must bear the following form: 

da \e+e~ -> JM + X] , , . o . 
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It is interesting to note that, after the suitable reduction, the squared matrix clement for 
this reaction can depend upon the polar angles only through rather limited combinations- 1, 
cos^ 6*, cos 6*1 cos 6*, and cos^ 6*1, respectively. Therefore, to arrive at the expression indicated 
in (42), suffices it to know the following integrations over d^^: 



L 



/•2n 



d^l = 2tt, 



d^l cos 9i — 27r cos 9 cos ©* , 



27r 



d^l COS 9l = 27r 



^ (1 + cos' ^) + - ^ ®i) (1-3 cos' 9) 



(43a) 
(43b) 
(43c) 



Thus we arc reassured that only the zeroth and second powers of cos 6 are allowed to appear 
in the double differential distribution for J/ip production, in conformity to (42). It may 
be also worth pointing out that, one great simplification can be made insofar as only the 
differential energy spectrum of unpolarizcd J/t/j is concerned. In this case, the second term 
inside the square bracket in (43c) can be discarded, since its contribution vanishes upon 
integration over 9. 



VI. INCLUSIVE J/ip DISTRIBUTIONS AT B FACTORY 

In this section, we report our results of order- relativistic correction to inclusive J/ip 
production associated with non-cc states at B factory at ^/s = 10.58 GeV. We investigate 
its impact on the integrated cross sections, and various types of distributions of J/ip at B 
factory, which is found to be modest. 

Very recently the order-v' correction to e+e~ — >■ J/i/^gg at B factory has also been studied 
by He, Fan and Chao [59], who have found similar magnitude of the relativistic correction to 
the integrated cross section for unpolarized J/ip. Since none of the differential distributions 
for energy, angular and polarization of J/ ip have been explicitly given in [59] , it is not possible 
at this stage to make a detailed comparison between our results and theirs. Nevertheless, 
these authors chose to expand Eg around rric in powers of q' in the amplitude, which is 
opposite to the strategy employed in this work. Most notably, it seems that relativistic 
correction effects associated with the three-body phase space has been neglected in [59], 
thus the exact agreement between our results and theirs will not be expected. 



A. Choice of input parameters 

To make concrete predictions, we need specify various input parameters, in particular 
the corresponding NRQCD production matrix elements. The LO NRQCD matrix element 
(Of^^), together with some specific combination of coupling constants and mass scales will 
be frequently encountered in many expressions for J/tjj production rate. For notational 
compactness, it is thus convenient to lump them into a single factor: 

= ^W^<f (44) 
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We take Mj/^ = 3.097 GeV, ^/s = 10.58 GeV at B factory energy, thus fix r = Mj^^/s = 
0.0857. For the NRQCD matrix elements, we quote the values extracted from the recent 
Cornell-potential- model-based analysis [54] : 

(O//'^) = 0.440^HI GeV^ (45a) 
{v')j/^ = 0.2251°;^°^. (45b) 

Note the uncertainties affihated with each NRQCD matrix elements are quite sizable. With 
resort to the G-K relation (18b), one finds that this value for {v'^)j/^ corresponds to the 

charm quark pole mass = 1.39 ± 0.06 GeV. 

Targeting at a better accuracy, we also including the running effect in the electromagnetic 
coupling, i.e., we take the fine structure constant to be a(\/i) = 1/130.9, rather than the 
commonly used 1/137 [23]. For the strong couphng constant, we take the central value ag 
equal to 0.21, corresponding to choosing the renormalization scale /j at v'i/2. The corre- 
sponding uncertainty is estimated by varying this coupling between 0.17 to 0.26, obtained 
by sliding the renormalization scale ii between ^/s and \/i/4 [23]. 

With all these parameters specified, we find 

ao - O.ldQt'oml Pb- (46) 
The attached error comes from the uncertainties of ag and of the NRQCD matrix element 



B. The integrated production rate for J ftp associated with Ught hadrons 

The lowest-order NRQCD predictions to the J/ip associated production rate have been 
available for a long while [25, 28-31]. For convenience of the reader, we collect in Appendix A 
the expressions for the energy distribution of (un)polarized J/t/j, at LO as well as at NLO 
in f ^ . 

One may attempt to directly integrate the differential J/ip spectrum over the entire J/ijj 
energy range to deduce the integrated J/ip cross section. Unfortunately, it seems rather 
difficult to obtain the analytic expression by this way, even for the LO cross section. Fortu- 
nately, to this purpose, it is much more advantageous to carry out the 3-body phase space 
integration in a route as specified in (23), where the corresponding integration boundaries 
become simpler. After some straightforward calculations, we find that the LO integrated 
cross section for the unpolarized J/ip can be put in the following compact form: 



f 2 — r — 12r^ -|- Sr^ 
+ Xiight] = (To I 2{l-r)^ arctanhVl - r 

4-9r + 8r2 , . 5 - 14r + 3r2 9(1 - 2r + 2r2) 

+ (i_,)3/2 arctanh^r^+ _ Inr- _ 



(47) 



We note that the analytic expression for (T*^°^ has already been available recently [44]. Our 
expression is in agreement with equation (2) in [44], but appears to be simpler 



This can be mainly attributed to the fact that a pair of dilogarithms appearing in their formula can 
actually be transformed away, by exploiting a sequence of identities about dilogarithms. 
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It is interesting to examine the asymptotic behavior of (47) in the high energy hmit 



1 f 1 \ 9 

- InV + In 2 In r + In^ 2 + 4 In 2 h 0(r InV) 

4 V 2 / 2 



(48) 



Beside the power-law scahng contained in uq (oc 1/s^), the asymptotic behavior of the total 
cross section is dominated by the double logarithm term. This expression is superficially 
analogous to the NLO perturbative correction to the exclusive double- charmonium produc- 
tion process e+e~ J /ip + ric, which also exhibits a double logarithm scaling [21]. 

It is also of some interest to examine the opposite limit r — )■ 1, in which the J/ijj is 
produced just above the kinematic threshold. The total cross section in this limit vanishes 
as (To(l — r)/3 + 0((1 — r)^), which may reflect that gluon radiation off the heavy quark is 
greatly damped in very restricted phase space. 

Substituting r = 0.0857 and the value of ctq given in (46) into (47), or equivalently, 
numerically integrating the spectrum (Al) over the entire J/ip energy, we find the LO 
prediction to the total cross section for J/ip associated with non-cc states at B factory is 

[ J/^ + Xught] = 0.200;°;^i pb. (49) 

The error is solely due to the uncertainty in (Tq. 

We now turn to the order-f ^ contribution to the integrated cross section for producing 
the unpolarized J/ip. Unlike in the LO case, the corresponding analytical expression is too 
complicated to be presented here, and we are content with providing numerical result only. 
Taking the value of (Jq from (46), together with the ratio of the NRQCD matrix elements 
{v'^)j/ip in (45b), integrating the order-v^ correction to the spectrum (A2) over the entire 
J/ip energy range, we find 

a(^) [j/ij + Xught] = o.oeilH'o pb, (50) 

where the attached error is due to the uncertainties in &q and in (f ^) j/^. It is clear to 
see that for the central value of the predictions, the inclusion of the order-i>^ correction 

enhances the LO result by about 30%. This seems in conformity to the naive expectation 
about the size of relativistic correction for charmonium system. It is interesting to note 
that, the central value of the relative order-w^ contribution seems even slightly larger than 
the recently-computed NLO perturbative correction, which enhances the LO result by about 
20% [43, 44]. But fairly speaking, the effects of both types of corrections are not significant. 
The sum of (49) and (50) turns to be 

(^(0) + + Xu,nt] = 0.2611^:?™ pb. (51) 

Compared with the latest BELLE measurements for prompt J/ip production rate associated 
with light hadrons, (lb), we find rough agreement between (51) and the data with large 
uncertainty. This agreement could be even more satisfactory if further including the NLO 
perturbative correction and the feeddown contribution from higher charmonium states. 

In light of this rough agreement achieved by the color-singlet contribution alone, one 
important question is to ask how much room is left for the color-octet contribution to 
inclusive J/ip production at B factory. It seems fair to state that earlier estimates of its 
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contribution [36, 37, 39] may turn out to be overly optimistic. Nevertheless, we would like 
to caution that, our predictions, both the LO one in (49), and the NLO one in (50), are 
subject to large theoretical uncertainty, so we are unable to draw any firm conclusion about 
the actual size of the color-octet contribution. 

C. Energy spectrum of unpolarized J/V' 




FIG. 2: The energy spectra of the unpolarized J/ijj (left panel) and longitudinally polarized J/'0 
(right panel) associated with light hadrons at the energy of B factory. The dot-dashed curve 
represents da^^^/dz, the dashed curve represents da^'^^/dz, and the solid curve represents their 
sum. For simplicity, in all the figures in this work, we have taken only the central values of the 
input parameters and not drawn the error band. 



Aside from the total production rate of J/ip, it is also useful to look closely into the 
differential observable. As a matter of fact, B factory experiments have already measured 
various types of J/ ip distributions. In this subsection we investigate the effect of first-order 
relativistic correction for the energy distribution of unpolarized J /ip. In Figure 2, we display 
the energy spectrum of unpolarized J/ip at i?-factory energy, including both the LO result 
and the first-order relativistic correction. 

As one can tell from Figure 2, the LO distribution admits a finite limit when the J/V' 
energy approaches its maximum: 



dz 



[e+e ^ J 1^1) + gg] 



2r 



CTo- 



(52) 



A novel feature of relative order-f ^ contribution is that, as can be clearly seen in Figure 2, 
the spectrum has a sharp rise near the very upper end of the J /if) spectrum. After some 
straightforward manipulation on the analytic expression of order-f ^ contribution, which is 
recorded in equation (A2), we find the following limiting value near the endpoint: 



da^ 

dz 



[e+e -> J/%1) + gg] 



>l+r 



X 



9 - 23r - lOr^ - 12r(l + r) In [ 
(1 -r)2 



r (1+r— z) " 
(l-r)2 . 



(53) 
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Clearly the endpoint singularity is of the form ln(l + r — 2;). 

The logarithmic divergence near the endpoint is not something new. It simply signals the 
breakdown of the NRQCD expansion near the kinematic boundary, as a result we should no 
longer trust our prediction in this region. Recall that for the NLO perturbative correction to 
the same process, the logarithmic singularity of hi(l + r — z) is also expected to appear near 
the maximum of J/ip energy [45, 46]. However, it is worth mentioning that, the ln(l + r — 
has rather different origin for both types of corrections. For the NLO perturbative correction, 
the ln(l + r — ^) term should be attributed to the coUinear singularity associated with the 
gluonic jet recoiling against J/t/j. The reason is that, at LO in v, the soft gluon cannot 
resolve the color-singlet cc pair (color-transparency), as a result the net contributions from 
soft gluons cancel, so the logarithm can be only of the collinear origin However, for the 
contribution from relativistic correction, this endpoint singularity comes from the region 
where one of the recoiling gluon becomes soft. Since we have gone beyond the LO in v, 
the color-singlet cc pair could still develop a nonzero color dipole, therefore it may strongly 
interact with the soft gluons. Therefore it is natural to identify this resulting ln(l + r — z) 
with the soft origin. It is interesting to ask whether the method presented in [45, 46], which 
combine NRQCD and the soft-coUinear effective theory, can be generalized to resum those 
types of logarithm in (53) to all orders in as, to render the J/ip energy spectrum well-behaved 
near the end point region. 

Note this endpoint singularity is integrable, therefore wc are still able to obtain a finite 
order-f^ correction to the integrated cross section (see (50)). This is similar to quarkonium 
semi-inclusive radiative decay J/ip — >■ + X, where the order-^;^ correction to the photon 
spectrum also develops an integrable endpoint singularity. Nevertheless in that 
relative order f^, the photon spectrum would develop a linear infrared divergence near the 
end point, which results in a logarithmic divergence for the integrated decay rate [57]. It 
is the color-octet mechanism that should be invoked to tame this infrared divergence. In 
our case, we expect the exactly same pattern will occur. That is, at 0{v^), the J/i/j energy 
spectrum would develop a linear endpoint singularity, consequently the integrated cross 
section would contain a logarithmic infrared divergence, which must in turn be cured by 
including the color-octet contribution. 



D. Polarization distribution of J/ij^ 

Babar and BELLE collaborations can also determine the polarization of J/ip as a func- 
tion of its energy by measuring the muons' angular distribution from J/xjj ^ l^'^l-i' - The 
commonly used polarization parameter is defined by 

. , da/dz — Sdai/dz , ^. 

'^{^) = , , , ; — n—, (54) 

^ ^ da/dz + daL/dz ^ ' 

where duj^jdz signifies the differential cross section for producing a longitudinally-polarized 
J li). a — 1 and —1 correspond to 100% transversely- and longitudinally-polarized, whereas 
a — corresponds to 100% unpolarized. 



It seems enlightening to contrast the single collinear logarithm associated with the color-singlet channel 
at LO in v with the Sudakov double logarithm associated with the color-octet channel [39] . 
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To deduce the function a{z), it is necessary to know the expression for dai/dz. The 
analytical expressions for this distribution, at both LO and NLO in v"^, have been given in 
Appendix A. Moreover, both the LO and NLO contributions to the energy spectrum for 
the longitudinally-polarized J/ip at factory is shown in Figure 2. 

Let us first investigate the integrated cross section for producing a longitudinally polarized 
J/ip. As in the unpolarized case discussed in Section VI B, if one carries out the 3-body 
phase-space integration following the order specified in (23), the LO integrated cross section 
for the longitudinally-polarized J/ip can also be put in a closed form: 

[J/V'IA = 0) + Xu,Ht] = ao { ^ ~ ~4(tl~ ^ arctanh^yr^ 

4 - 4r + - 3r^ 2 - lOr 7r^ - 6r^ + 3r^ , 

vU^ arctanh VI - r H — In r 

2(i-r)3/2 4(1 -r)2 

-I- ^ I Li2(v^) — Li2(— v^) ~ arctanhvr mr j 

+ — 2-12V^+3r + 2r2 + 3r2-6r2 + — ^ }, 55 

16 \ / 4(1 — r) ) 

where Li2 stands for the dilogarithm. This analytic expression has not been known previ- 
ously. A nontrivial check of the correctness of this result is to examine its threshold behavior. 
In the limit r — )■ 1, crj^^ approaches zero as (Tq^ + 0{{1 — r)^), as expected, one-third of 
that for the polarization-summed case. 

It is of interest to ascertain the asymptotic behavior of a^^^ in the high energy limit 
> Mj/^: 

[J/V'(A = 0)+XH,ftt] 



^ InV + Q - In 2^ In r + In^ 2 - 2 In 2 + 2 + y + 0{^/^) 



(56) 



lni + 61n2-f 



Note the leading double logarithm appearing here has the same coefficient as in the 
polarization-summed case, (48) One can then readily infer the asymptotic behavior of the 

transversely-polarized J/xjj production rate: a!^^ = — cr^^ — >■ <^o 
only exhibiting a single-logarithm scaling. 

Substituting r = 0.0857 into (55) and using the value of &q given in (46), or straightfor- 
wardly integrating the spectrum (A3) over the entire J/ip energy numerically, we find the LO 
prediction to the integrated rate for producing longitudinally-polarized J/ip in association 
with non-cc states at B factory to be 

[J/V' + Xu,Ht] = 0.128t°;°i pb. (57) 

The error originates solely from the uncertainty in (Jq. According to (54), and using the 
central value of (49) and (57), we find the a — —0.56 averaged over the entire J/-^ energy 
range. 



In contrast to (48), the leading correction to this asymptotic expression is of relative order instead 
of 1/s. 



26 



For the order-?;^ contribution to ai, the corresponding analytic expression is too involved, 
if not impossible, to deduce, so we are content with providing numerical result only. Using 
as given in (45b), integrating the order-v^ correction (A4) over the entire J/ijj energy 
range, we get 

(jf [J/iP + Xiight] = 0.037tH24 Pb. (58) 

The attached error comes from the uncertainties in (Jq and {v'^)j/^. For the central values of 
the predictions, inclusion of the order-f ^ correction enhances the LO cross section by about 
29%, which has a very similar magnitude of enhancement as for the unpolarized J/ijj. This 
is again in accordance with the naive expectation about the size of relativistic correction for 
charmonium system. 

The sum of (57) and (58) is 

(4°^ + '^?^)[J/^ + XugHt] = 0.165t'o^l pb. (59) 

Including the order-f^ correction, the central value of the average polarization variable a 
shifts from —0.56 to —0.55. Hence the relativistic correction has a rather minor effect in 
changing the polarization of J/ip. 

Now let us examine the differential distribution daf^ /dz. As can be seen from Figure 2, 
or can be directly inferred from (A3), the LO distribution has a finite limit when the energy 
of the longitudinally-polarized J/ if) approaches its maximum: 



da 



(0) 
L 



dz 



[e+e-^ J/V^(A = 0)+55] 



ao-^. (60) 



z— >-l+r 



From (52) and (60), it is ready to see that, at the endpoint 2; = 1 + r, the LO polarization 
variable, a^^\ approaches the constant —j^- 

As can be seen in Fig. 2, the order-v^ correction to the energy spectrum of the lon- 
gitudinally polarized J/ijj also diverges logarithmically near the upper end. After some 
manipulation on equation (A4), we find the following limiting behavior: 



da 



(2) 



dz 



[e+e-^ J/V^(A = 0)+^^] 



^0 {v'] 



3 



5- 17r + 4r2-8rln [^^±^1 
X (1-.). ■ (^1) 

The situation very much resembles that for the unpolarized J/V'- One can refer to the 
paragraphs after (53) for similar discussions. 

In Figure 3, we also display how the polarization parameter a varies with the J/tl) energy. 
Clearly, the inclusion of relativistic correction seems to have a minor impact in most of the 
region of except increasing it modestly near the upper end. 



E. Angular-Energy distribution of J/i/' 

Experimentally it is also possible to measure the production rate for J /ip in e^e~ annihi- 
lation that is differential in cos^, the cosine of the angle between the momentum of J/%1) and 
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a (z) 



-0.6 




FIG. 3: Profile of tlie J /ip polarization parameter a{z) in associated production with light hadrons 
at -^s = 10.58 GeV. The dot-dashed curve represents the leading order prediction q^''-', whereas 
the solid curve represents the corresponding one including the 0{v'^) effect. 



the incident e~ beam in the laboratory frame. It is thus theoretically interesting to study the 
differential angular distribution of J As pointed out in (42), for inclusive J /ip production 
in e^e~ annihilation, general consideration constrains the double differential distribution of 
the following form [28]: 



dzd cos 9 



[e+e- ^ J ftp + X] = S^'^iz) [1 + A^'^iz) cos' 



(62) 



where A{z) is a angular parameter that satisfies |A(-2)| < 1. 

The analytic expressions at LO in v, S^'^^z) and A^'^\z) have been known long ago. The 
closed forms of the order-f' contributions, S^'^\z) and A^'^\z), are derived in this work for 
the first time. For completeness, we reproduce all of them in Appendix A. From Eqs. (AG) 
and (A7), one finds that the LO double differential spectrum admits a finite limit near the 
upper end [36]: 



dz d cos 9 



[e+e ^ J/i/j + gg] 



3 (To / 1 + r 



>l+r 



cos^9 



(63) 



which implies that at the endpoint z = 1 + r, A^'^^ = — j:^. 

In Figure 4, we display the angular function A{z) at energy of the B factory, ^/s = 10.58 
GeV. Both the LO prediction and that including the first-order relativistic correction are 
shown. Note that the correct A{z) incorporating the order-f' effect is given by [28] 



^(0)(;,)A(0)(;,)+g(2)(^)^(2)(^) 

S(^\z) + S(^'){z) 



(64) 



As can be seen in Figure 4, including the first-order relativistic correction seems to have 
modest effect, which only slightly softens the angular distribution near the upper end. 

Next let us inspect the end-point behavior of the functions S^'^\z) and A^'^\z). After some 
straightforward algebra from Eqs. (A8) and (A9), we find the following limiting behaviors 
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A(z) 




FIG. 4: Profile of the J /ip angular distribution parameter A{z) in associated production with light 
hadrons at = 10.58 GeV. The dot-dashed curve represents the leading-order prediction A^^\ 
whereas the solid curve represents (z) defined in (64), which has included the order- 1;-^ effect. 



when z approaches its maximum: 



a, {v^)j/.^ 5 - 16r - 5r' - 2r(5 + 3r) In [^ii^] 

'1-r) 



A / 1 W 



4 1-r 

These emerging logarithmic divergences near the endpoint simply reflects that the differential 
cross section diverges in that region. However, according to Eq. (64), the angular distribution 
(z), defined ratio, still remains a finite and smooth function near the very upper 
end. 

It is worth mentioning that BELLE collaboration has recently measured the average 
angular variable for the J/ip production in association with noncharmful states, Aexp = 
5.2l2i(0.3) [40]. Theoretically, it is straightforward to define the corresponding A by inte- 
grating (62) over z: 

"^""^^^ ■ [e+e- ^ J/^ + X,,,,,] = S^'^ [1 + A^'^ cos' 9] , (66) 



dcos9 



From Eqs. (A6) and (A7), and inserting r = 0.0857, we find the LO NRQCD prediction is 

A = —0.037. Notwithstanding the large experimental uncertainty, this prediction is in 
apparent disagreement with the Belle measurement, even the sign is opposite. Subsequent 
studies reveal that including the NLO perturbative correction does not help to resolve this 
discrepancy [42]. 

One may naturally wonder whether implementing the relativistic correction will bring 
the NRQCD prediction closer to the data or not. In analogy with (64), we introduce a new 
average angular variable that incorporates the 0(f ^) effect: 
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Starting from Eqs. (A8) and (A9), and adopting the central value of tabulated in 

(45b), we then find A = 0.0011, which now has the same sign as the measured value, though 
still differs considerably in the absolute magnitude. Therefore, including NLO rclativistic 
correction (plus perturbative correction) seems not to be sufficient to explain the data. It 
remains to be a challenge how to correctly account for the measured angular distribution in 
the J/ip + Xiigfit channel within the NRQCD factorization framework. 

VII. DISCUSSION AND SUMMARY 

In this work, we have introduced a somewhat heterodox NRQCD matching strategy, 
which is particularly suitable for calculating the rclativistic correction to (inclusive) quarko- 
nium production and decay processes with involved kinematics in the color-singlet channel. 
The great advantage of our approach over the orthodox matching strategy is that, it can 
take into account the rclativistic correction effect in the phase space integration with much 
ease, thanks to the Gremm-Kapustin relation. As a nontrivial application of this method, 
we have systematically investigated the relative order-^;^ correction to the inclusive J/ip pro- 
duction associated with hght hadrons at B factories. We have found that it can modestly 
enhance the lowest-order NRQCD prediction for the integrated J/ip cross section, about 
30% if we choose the rclativistic correction matrix element as specified in [54]. We find its 
impact on the J/ijj polarization and angular distributions is quite minor. The magnitude 
of the order-v^ correction seems to be comparable with that of the respective NLO pertur- 
bative correction. We would like to caution that, our predictions of the order-f ^ correction 
are likely subject to large theoretical uncertainty. In particular, some intrinsic uncertainty 
related to the relative order-v^ NRQCD matrix element seems to restrict our ability to make 
precise predictions for the rclativistic correction. Since the corrections computed in this 
work has only a modest effect, we feel unable to draw any sharp conclusion, especially for 
the actual size of the color-octet contribution. 

Of the special theoretical interest, is the logarithmic divergence near the upper end of 
the J/ijj spectrum found in this work for the order-i;^ contribution. It is desirable to extend 
the theoretical framework developed in [39, 45, 46] beyond the LO in i>, to see whether such 
type of soft endpoint logarithms can be resummed to all orders in ct^, to render the J/ijj 
energy spectrum well-behaved in the end point region. 

Another interesting direction is to incorporate the order-v*^ correction to the process 
considered in this work. We expect that the perturbative matching approach described 
in this work, after some straightforward extension, is well-suited to achieve this goal. At 
O(f^), the J/ip energy spectrum is expected to exhibit a linear divergence near the upper 
end point, and consequently, the integrated cross section will be logarithmically divergent. It 
will be interesting to see how the color-octet contribution from the ^Pj^^ NRQCD production 
operator is explicitly put into work to tame this infrared divergence. 
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Appendix A: Miscellaneous formulas for inclusive J/ip production associated with 
light hadrons in e+e" annihilation 

In this section, we collect the analytic expressions for various types of distributions for J/ip 
associated production with light hadrons in e^e~ annihilation. Each type of differential cross 
section is understood to contain two parts: da = da'^'^^ + da'^'^\ which represent the leading 
order contribution, and the contribution of relative order-v^, respectively. We emphasize 
that it is the physical J/ip mass, rather than the charm quark mass, that enters into the 
formulas of each part. 



1. Energy distribution of unpolarized J/tp 

The energy spectrum of unpolarized J/ip at LO in v reads: 

[e+e ^ J/iP + gg] 
1 

= (To 



(2-z)2(2-2r)3 

X ^{z - 2r)\/ z"^ - Ar [4(1 + 5r + 7r^ + 4r^) 

- 12(l + r)(l + 2r)2 + (13 + 14r)2^-4/] 

+ 4(1 + r - [2r(l - r)(l + 8r + 4r^) - 2r(5 - 2r - 6r^)^ 

. (l..-5.V].n( ;:^;:gIg }, (Al) 

where the quantity (Jq has been defined in (44). This expression agrees with the result given 
in [25], but differs from Ref. [28-30] by an overall constant. 
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The first-order relativistic correction to the energy spectrum of unpolarized J/tl) reads: 
^— [e+e ^ J/il) + gg] 

_ O'O (''^).//r 1 

~ 3 (2 - zY{z - 2rf 

X |(z - 2r)V-z2 - 4r [64r(3 + llr - 2r^ - 4r^ - 20r^ - 15r^) 

- 32r(22 + 30r - 21r^ - 91r^ - 89/ - 7r^)z 

- 16(1 - 48r - 9r^ + 171r=^ + 213r^ + 35r^)z^ 
+ 16(4 - 22r + 66r2 + 133r^ + 35r^)z^ 

- 4(18 + 15r + 170r^ + 70r^)z^ 

+ 4(ll + 15r + 16r2)z^-(ll-2r)^^] 

+ 4 [32r2(l + r)(3 + 9r - Qr'^ + 9r^ + 14r^ + 15r^) 

- 16r^(24 + 50r + 21r^ + 124r^ + 228r^ + 126r^ + 7r^)z 

- 8r(3 - 63r - 103r^ - 278r^ - 697r^ - 469r^ - 49r^)z^ 
+ 8r(9 - 50r - 186r^ - 549r^ - 477r^ - 77r^)z^ 

+ 2(2 - 37r + 248r2 + 993r^ + 1122r^ + 272r^)z^ 

- 4(2 - 2r + 124r^ + 192r3 + 71r^)^^ 
+ (7 + 26r + 140r^ + 87r^)z^ 

- (3 + 2. + 14.V].n( -^;:gI| )}. (A2) 
where has been introduced in (5). 

2. Energy distribution of longitudinally-polarized J/ip 

The energy spectrum of the longitudinally-polarized J /tl) at LO in v reads: 

^ [e+e- ^ JMA = 0) + ^5] 

^ "° (2 - zY{z - 2r)3(.^ - 4r) 1^" " 2r)V-^^ 
X [-8r^(9 + 9r + + 3r^) + 16r(2 + 8r + + ?,r^)z 

- 4(1 + lOr + 3r2 + Sr^)^^ + 4(1 - r + 2r2)z^ + ^^j 

+ 4(1 + r - 2) [-4r^(l - r)(9 + 2r + 3r^) + 8r^(2 - 4r - 3r^)z 

- 2r(l - 16r - 2r^ - 9r^)z^ - 2r(5 + 3r + 3r^)/ 

+ (l + '- + '-V]ln( ^_^^_^^,_J |. (A3) 
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The order-v^ correction to the energy spectrum of the longitudinally-polarized J /tl) is 



dz 



3 (2 - z)'^{z - 2r)5(z2 - 4r) 
X |(z - 2r)Vz^ - \r [-128r^(15 + 4r + 5r^ + - 4r^ - 15r^) 

+ 64r2(6 + 62r + SOr^ + 39r^ + 21/ - lOlr^ - 9r^);2 

- 32r2(6 + 122r + \h'ir^ + 215r^ - 247r^ - hlr')z^ 

- 16r(2 + 6r - Tlhr^ - 57lr=^ + 233r^ + \Mr')z^ 

- 16(1 + 7r + 51r^ + 358r^ + 16r^ - 83r^)^^ 
+ 8(6 + 45r + 205r^ + 113r^ - 41r^)^^ 

- 4(14 + 91r + 76r^ + 3r^)z^ 

+ 4(8 + 15r + 5r^)/ - (3 + \r)z^\ 

+ 4 [-64r^(15 + 9r + 3r2 - 2r^ - + ^r' + 15r^) 

+ 32r^(6 + 74r + 32r2 - 5r^ - 28r^ + 26r^ + 134r^ + 9r^)z 

- 16r^(108 - 9r - 166r^ - 118r=^ + 478r^ + 75r^)z^ 

- 8r2(6 + 116r + 71r2 + 492r^ + 784r^ - 840r^ - 261r^)^3 

- 8r(l - 21r - 195r^ - 437r^ - 947r^ + 308r^ + 2A3r^)z^ 

- 4r(2 + 63r + 497r2 + 1253r^ + 133r^ - 252r^);2^ 
+ 2(2 + 21r + 189r2 + 964r3 + 447r^ - 127r^)^^ 

- 4(2 + 16r + 91r2 + 86r^ - 
+ (7 + 46r + 48r2 + llr^)^^ 

3. Angular- energy distribution for unpolarized J/ip 

The doubly differential angular-energy distribution of unpolarized J/ip can be parame- 
terized in the following form: 

drr^^ 

[e+e- J/i/j + gg] = S^"^ {z) [l + A^"^ {z) cos' 9] , (A5) 



dzd cos 9 



where the superscript i = 0,2 represents the leading-order and first-order contributions in 
relativistic expansion. 
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The corresponding functions at LO in v read: 

3 (To 1 



4 (2 - - 2r)3(^2 _ 

X |(^ - 2r)Vz^ - 4r [ -4r(l + r)(3 + 12r + 13r^) 

+ 32r(l + r)(l + 3r)z + 4(1 - 7r - 12r2 + 2r^)z^ 

- 4(2 + r + 3r^)z^ + 7(1 + r)z^ -2z^] 

- 2(1 + r - [4r2(l - r)(3 + 24r + ISr^) - 8r\7 - 3r - 12r^);2 
+ 2 r(l - lOr - 27r^ + Ar^)z^ + 2r(7 + 7r - 6r^)^^ 

- (1-0(1 + 5.).'] ln(-^-^_j|. 



and 



S^'\z)A^'>\z) 



3 (To 



4 (2 - z)^{z - 2r)3(^2 _ 4r) 
X |(z - 2r)V-z2 -4r [4r(l + 5r + 19r^ + 7r^) 

- 96r2(l + r)z - 4(1 - 5r - 22r2 - 2r^)z'^ 

- 4r(7 + 3r)z^ + (5 + 7r)z^-2z^] 

+ 2(1 + r - [4r^(l + 7r)(l - r^) - 8r^(l + 3r)(l - 4r)z 

- 2r(l + lOr + 57r^ + Ar^)z'^ + 2r(l + 29r + Qr'^)z^ 

Note these expressions are exactly twice smaller than Eqs. (Ala) and (Alb) in Ref. 
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At the relative order v'^, the corresponding functions S{z) and A{z) are 



X <{z- 2r)Vz^ -Ar [ -64r2(9 + 35r - Sr^ - Ur^ - 57r^ - 45r^) 



+ 32r2(64 + 84r - 89r^ - 269r3 - 275r^ - 19r^)^ 
+ 16r(13 - 69r + Glr^ + 557r^ + 638r^ + 64r^)^2 

- 8r(110 + 161r + 559r^ + 725r^ + 3r^ - Ur^)z^ 
+ 4r(287 + 451r + 383r2 - 251r3 - 70r^)z^ 

+ 8(2 - 87r - 29r2 + 109r^ + 35/)^^ 

- 4(6 - 32r + 75r^ + 35r=^)/ 

+ 2{9 + 9r + 16r^)z'^-{5-r)z^] 

+ 2 [ -6Ar^{9 + 38r + 13r^ + ISr^ + 53r^ + BOr^ + 45r^) 

+ 32r^(70 + 166r + 123r2 + 340r=^ + 620r^ + 390r^ + 19r^)z 

+ 16r^(19 - UOr - 430r^ - 922r^ - 1927r^ - 1410r^ - 102r^)z^ 

- 8r^(124 - 219r - 1406r^ - 3154r^ - 2596/ - 139r^ + Ur^)z^ 

- 4r(18 - 193r + 840r^ + 2976r^ + 2400r^ - 327r^ - 98r^)/ 
+ 8r(29 + 28r + 355r2 + ISlr^ - 352r^ - 77r^)z^ 

+ 2(2 - Ulr - 131r^ + 293r^ + 1041r^ + 272r^)z^ 

- 4(2 - 34r + 71r^ + 196r=' + 71r^)z^ 
+ (7 - 5r + 147r^ + 87r^)^^ 



S(^\z) 



^0 

4 



(2 - z)^{z - 2r)5(^2 _ 



1 



(3 + r + Ur^)z^] In 



( 



z -2r + y/z"^ - 4r 
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and 



4 (2 - zY{z - 2rf{z'^ - 4r) 
X |(z - 2r)Vz'^ - 4r [64r2(3 + 17r - 8r^ - lOr^ - ll/ - 15r^) 

- 32r2(16 + 12r - 99r^ - 79r^ - 113r^ - r^)z 

- 16r(7 + 89r + 271r=^ + 335r3 + 370r^ + 32r^);2=^ 
+ 8r(90 + 595r + 789r^ + 775r=* + 161/ + 14r^)^^ 

- 4(1 + r)(8 + 325r + 836r^ + 32 Ir^ + 70r^)z^ 
+ 8(10 + 129r + 291r2 + 141r^ + 35r^)z^ 

- 4(18 + 104r + 119r^ + 35r^)/ 

+ 2(l + r)(17+16r)/- (7-r)^^] 

+ 2 [64r^(3 + 18r + 15r^ + 30r^ - 25r^ + 8r^ + 15r^) 

- 32r^(18 + 98r + 201r2 + 28r=^ + 36r^ + 162r^ + r^)z 

- 16r^(9 - 12r - 490r^ - 566r^ - 389/ - 710/ - 34/)^^ 

+ 8/(36 - 257r - 1410/ - 1670/ - 1980/ - 193/ - 14/)/ 
+ 4r(14 + 221r + 1360/ + 3208/ + 3800/ + 595r^ + 98/)/ 

- 8r(35 + 300r + 817/ + 1203/ + 284/ + 77/)/ 
+ 2(2 + 219r + 1177/ + 1973/ + 669/ + 272/)/ 

- 4(2 + 82r + 275r2 + 124r^ + 71r^)/ 
+ (7+119r + 119/ + 87/)/ 



(3 + 5r + 14/)/] In 



-2r + ^/z 



4r 



(A9) 



— 2r — v-z 4r^ 

Integrating Eq. (A5) over the polar angle 9 from to tt, we arrive at the following identity: 



dz 



[e+e- ^ J/ip + gg] ^2S^\z) 



1 + -A«(z) 



(AlO) 



where da'^^^dz represent the energy distributions for unpolarized J/4', which have been 
given in Eqs. (Al) and (A2). This relation can serve as a consistency check of our results. 
We have explicitly verified that our expressions obey this relation for both i = and 2. 



Appendix B: Equivalence between our matching method and the "orthodox" one 

It is curious to ask, whether the inclusive J/ip production rate derived from our matching 
procedure, can be translated into a more orthodox form, that is, everything is expressed in 
terms of charm quark mass rather than the charmonium mass. That corresponds to what 
would be obtained from a literal matching method. As we shall see, this is possible only 
for the integrated cross section. And we like to stress, there should be no any theoretical 
ambiguity and confusion for the relativistic correction contribution at this level. 
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To better orientate ourselves, let us begin with a one-dimensional toy integral: 



dtW{t,y) = / dtW{t,0) + y\ / dtW;^{t,Q) 

g{x,y) Jg(x,0) yjg{xfl) 

+ W{f{x,Q),Q) fy{x,Q) - W{g{x,Q),Q) g'y{x,Q)] + 0{y^), (Bl) 

where we have assumed the integrand W is regular at the end points of the integral and used 
the shorthand fy{x,0) = df{x,y)/dy\y=Q. The final result of the integral in the left-hand 
side will be a function of x and y. Here y is assumed to be a small constant, and it is 
assumed that both the integrand W and the integration boundaries /, g depend on y. Our 
goal is to reexpress the original integral in a Taylor-series in y. Since in many situations, 
the closed form for such integral is presumably not available, or at least difficult to obtain, 
it is thus desirable to find a general numerical recipe to accomplish this expansion. 

In the right-hand side of (Bl), we give the intended answer for this expansion through 
the linear order in y. The leading term is obtained by neglecting y simultaneously in the 
integrand and integration boundaries. The coefficients of order y come from either expanding 
the integrand or taking into account the correction to the integration boundaries. 

The goal is to reexpress our "leading-order" cross sections in terms of a new series in- 
cluding the first-order relativistic correction, with all the occurrences of Mj/^ replaced by 
2mc in a consistent way. Note both the matrix element squared and the boundaries of the 
phase space integral depend on v'^ implicitly through Mj/^. Clearly, is the counterpart 
of y in (Bl) that acts as the small expansion parameter. One can utilize (Bl) to work out 
the desired expanded form. 

For concreteness, we take the unpolarized J/ip energy distribution as an example. The 
LO energy spectrum of J/ip derived from our matching method has been given in Eq. (Al): 



dz 

2567r(ecQ;Q;<j)^ 



27M 



[e+e ^ J/tP + gg] 



1 



X 



I (2; - 2r)Vz'^ - 4r 



(2 - z)^{z - 2r)3 
4(1 + 5r + 7r^ + 4r^) 



- 12(1 + r)(l + 2r)z + (13 + Ur)z^ - iz^ 



+ A{l + r-z) 
+ (1 + r -5r'^)z'^ 



2r(l - r)(l + 8r + Ar^) - 2r(5 - 2r - 6r^)z 



In 



-2r + \/z^ - 4r 
— 2r — \J z^ — At 



(B2) 



For clarity, here we abandon the use of the abbreviation uq and supply the complete expres- 
sion of the prefactor. 

In accordance with (Bl), we may reexpress the integrated cross section of (B2) as a sum 
of the following three terms, each of which now depends on the charm quark mass rather 
than the J jil) mass: 



/ ^^^r~ = / / ~i — + (7(2''^ +0(vV). 

dz dz dz 



(B3) 
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where tq = Upon expanding (B2), we need replacing every occurrence of Mj/^p 
the combination of rric and (t'^)j/.!/) through the G-K relation (18b): 



with 



M 



ro[l + {v')j/^ + 0{v')], 



(B4a) 
(B4b) 



In the resulting new expression, we only need retain those terms at most of order v'^. 

The first term in the right side of (B3) constitutes the leading contribution, the second 
one comes from the expansion of the integrand, and the third one arises from the correction 
due to integration boundaries. Their explicit expressions are 



dz 



dz 

647r(ecaQ;s)^ 



dz 



27m,s2 ^ '(2-2)2(^-2r)V^2_4r 
X {{z-2r)[ -32r\5 + 17r - 4r^ - 12r^) 

- 16r(l - 17r - Sr^ + 32r^)z 
+ 8r(15 + 7r + 17r^ - \2r')z^ 

- 4(3 + 47r + llr^ - 32r^)z^ 

+ 2(10 + 51r - Mr^)^^ - 13(1 + 2r)z^ + \z^\ 

+ 4V^2 - 4r [4r^(5 + 24r + 3r^ + 8r^ + 12/) 
+ 2r(l - 36r - 45r^ - 32r^ - 56r^)^ 
+ 2r(l + 30r + + 53r^)z2 

- (1 + 2r + 34r2 + 55r^)^=^ 

— 2r + V-z^ — 4r 



(B5a) 
(B5b) 



+ 



(1 - r + ISr^)/] In 
5127r(ecQ; Q;s)^m, 



27s3 



J/V 



z-2r- 
l + 2r 



4r 



r—>-ro 



(B5c) 



where {v(^^) is given in (4b). Needless to say, the new LO term is exactly of the same 
functional form as the old one in (B2), except Mj/^ everywhere replaced by 2mc. For the 
newly generated relativistic correction pieces cr^^o) ^^d 5-^^''^ , one does not need to carefully 
distinguish tq and r in them, since the induced error would be of order which is beyond 
the intended accuracy of this work. 

All these three terms, in combination with (A2), the genuine 0{v^) contribution in our 
matching approach constitute an alternative but equally valid prediction to the integrated 
J/i/j cross section that is accurate at relative order v'^. Since the expression for the integrated 
J/ip production rate, when everything is expressed in term of rric, has no any ambiguity 



Note we can carelessly replace Mj/^ by 2mc in (A2), and, in the corresponding phase-space integral 
boundaries, since the induced error would be O(u^). 
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through 0{v'^), it can be used to check the correctness of the calculation performed in an 
"orthodox" matching method {e.g., see [59]). 

To clearly see how (B3) works, we can take advantage of our analytic knowledge for the 
integrated J/ip cross section at 0{v'^). Directly Taylor expanding (47) around r = Tq to first 
order in r — Tq, we find 

~(2a) , ~i2b) ^ l287T{e,aa,y j,^ f 2 - 9r + SQr^ - 2Sr' + 8r^ 
^'^ 27mes2 ^'^ ^\ 4(1 -r)3 

, 2 n 3 - 4r - + 4r^ , 

X arctanh yl — r -\ ; arctanhy 1 — r 

(1 - rp^ 

5-llr + 33r2-3r3 11 - 19r - 46r2 + ISr^ 

+ — 4(r^7p — 4(1^7)2 /■ 

We have numerically compared (B6) with the sum of (B5b) and (B5c) upon integration over 
the full range of z, and indeed found exact agreement. 

We have also numerically checked that, both sides of (B3), assuringly, do agree with each 
other at the integrated level, up to an error of order f ^ 

Finally, it might be worth mentioning that, the differential distribution da^'^^'^dz diverges 
at both upper and lower ends of z (albeit being the integrable singularities): 



dz 
da^^''^ 



..2V. 27m..2 \l-Vr)V.^-4r' ^^'"^ 



dz 



\2 



2 



(1 - 3r)(l + llr + 6r2) + 8r(l - 3r - r^) In r^^+V^ I 

X LJlzlLJ., (B7b) 

(l-r)3 ^ ' 

These artificial end-point singularities affiliated with the relativistic correction to J/ip energy 
distributions, especially the one appearing at the lower end, are clearly at odds with one's 
expectation and certainly not favored by the data. This may signal that, even if feasible, 
it is not of much benefit to perform the NRQCD matching in a strictly orthodox ansatz. 
Instead the matching method described in this work seems much more satisfactory. 



[1] For a comprehensive, but a slightly outdated review on quarkonium production, see N. Bram- 

billa et al. [Quarkonium Working Group], arXiv:hep-ph/0412158. 
[2] For a latest review, see G. T. Bodwin, Quarkonium Production and Decay: NRQCD Confronts 

Experiment, talk given at the KITPC-EFT-2009 program. The content can be downloaded at 

the following URL: http://www.kitpc.ac.cn/program.jsp?id=PE20090720&:i=sched. 
[3] G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51, 1125 (1995) [Erratum-ibid. D 

55, 5853 (1997)] [arXiv:hep-ph/9407339]. 



For J/tjj production at B factory, the relativistic correction stemming from expanding the phase space 
boundaries, a-^^^\ makes negligible contribution due to the additional suppression by m^/s. 



39 



[4] E. Braaten and S. Fleming, Phys. Rev. Lett. 74, 3327 (1995) [arXiv:hep-ph/9411365]. 
[5] P. L. Cho and M. B. Wise, Phys. Lett. B 346, 129 (1995) [arXiv:hcp-ph/9411303]. 
[6] J. M. Campbell, F. Maltoni and F. Tramontano, Phys. Rev. Lett. 98, 252002 (2007) 
[arXiv:hep-ph/0703113]. 

[7] P. Artoisenet, J. P. Lansberg and F. Maltoni, Phys. Lett. B 653, 60 (2007) [arXiv:hep- 
ph/0703129]. 

[8] B. Gong and J. X. Wang, Phys. Rev. Lett. 100, 232001 (2008) [arXiv:0802.3727 [hep-ph]]. 
[9] B. Gong and J. X. Wang, Phys. Rev. D 78, 074011 (2008) [arXiv:0805.2469 [hep-ph]]. 
[10] B. Gong, X. Q. Li and J. X. Wang, Phys. Lett. B 673, 197 (2009) [arXiv:0805.4751 [hep-ph]]. 
[11] G. C. Nayak, J. W. Qiu and G. Sterman, Phys. Lett. B 613, 45 (2005) [arXiv:hep-ph/0501235]; 
Phys. Rev. D 72, 114012 (2005) [arXiv:hep-ph/0509021]; Phys. Rev. D 74, 074007 (2006) 
[arXiv:hep-ph/0608066]. 

[12] G. C. Nayak, J. W. Qiu and G. Sterman, Phys. Rev. Lett. 99, 212001 (2007) [arXiv:0707.2973 

[hep-ph]]; Phys. Rev. D 77, 034022 (2008) [arXiv:0711.3476 [hep-ph]]. 
[13] K. Y. Liu, J. P. Ma and X. G. Wu, Phys. Lett. B 645, 180 (2007) [arXiv:hep-ph/0601215]. 

[14] X. G. Wu and Z. Y. Fang, arXiv:0904.3206 [hcp-ph]. 

[15] M. Beneke, I. Z. Rothstein, and M. B. Wise, Phys. Lett. B 408, 373 (1997) [arXivrhep- 
ph/9705286]. 

[16] T. Mannel and G. A. Schuler, Z. Phys. C 67, 159 (1995) [arXiv:hep-ph/9410333]; 
T. Mannel and S. Wolf, arXiv:hep-ph/9701324; 

L Z. Rothstein and M. B. Wise, Phys. Lett. B 402, 346 (1997) [arXiv:hep-ph/9701404]; 

M. Bcncke, G. A. Schuler and S. Wolf, Phys. Rev. D 62, 034004 (2000) [arXiv:hep- 

ph/0001062]. 

[17] Z. B. Kang, J.-W. Qiu and G. Sterman, in preparation. The preliminary 
content can be found in the talk delievered by J.-W. Qiu at the KITPC- 
EFT-2009 program, PQCD factorization for heavy quarkonium production, 
http:/ /www.kitpc.ac.cn/program.jsp?id=PE20090720&i=sched. 

[18] K. Abe et al. [Belle Collaboration], Phys. Rev. Lett. 89, 142001 (2002) [arXiv:hep- 
ex/0205104]. 

[19] E. Braaten and J. Lee, Phys. Rev. D 67, 054007 (2003) [Erratum-ibid. D 72, 099901 (2005)] 
[arXiv:hep-ph/0211085]; 

K. Y. Liu, Z. G. He and K. T. Chao, Phys. Lett. B 557, 45 (2003) [arXiv:hep-ph/0211181]. 
[20] Y. J. Zhang, Y. J. Gao and K. T. Chao, Phys. Rev. Lett. 96, 092001 (2006) [arXiv:hep- 
ph/0506076]. 

[21] B. Gong and J. X. Wang, Phys. Rev. D 77, 054028 (2008) [arXiv:0712.4220 [hep-ph]]. 
[22] Z. G. He, Y. Fan and K. T. Chao, Phys. Rev. D 75, 074011 (2007) [arXiv:hep-ph/0702239]. 
[23] G. T. Bodwin, J. Lee and C. Yu, Phys. Rev. D 77, 094018 (2008) [arXiv:0710.0995 [hep-ph]]. 
[24] T. V. Uglov, Eur. Phys. J. C 33, S235 (2004). 
[25] W. Y. Keung, Phys. Rev. D 23, 2072 (1981). 

[26] J. H. Kuhn and H. Schneider, Z. Phys. C 11, 263 (1981); Phys. Rev. D 24, 2996 (1981). 
[27] V. M. Driesen, J. H. Kuhn and E. Mirkes, Phys. Rev. D 49, 3197 (1994). 

[28] P. L. Cho and A. K. Leibovich, Phys. Rev. D 54, 6690 (1996) [arXiv:hep-ph/9606229]. 
[29] F. Yuan, C. F. Qiao and K. T. Chao, Phys. Rev. D 56, 321 (1997) [arXiv:hep-ph/9703438]. 
[30] S. Baek, P. Ko, J. Lee and H. S. Song, J. Korean Phys. Soc. 33, 97 (1998) [arXivihep- 
ph/9804455]. 

[31] K. Hagiwara, E. Kou, Z. H. Lin, C. F. Qiao and G. H. Zhu, Phys. Rev. D 70, 034013 (2004) 



40 



[arXiv:hcp-ph/0401246]. 

[32] A. V. Bcrczhnoy and A. K. Likhoded, Phys. Atom. Nucl. 67, 757 (2004) [Yad. Fiz. 67, 778 

(2004)] [arXiv:hep-ph/0303145]. 
[33] D. Kang, J. W. Lee, J. Lee, T. Kim and P. Ko, Phys. Rev. D 71, 094019 (2005) [arXiv:hep- 

ph/0412381]. 

[34] B. Aubert et al. [BABAR Collaboration], Phys. Rev. Lett. 87, 162002 (2001) [arXiv:hep- 

ex/0106044]. 

[35] K. Abe et al. [BELLE Collaboration], Phys. Rev. Lett. 88, 052001 (2002) [arXiv:hep- 
ex/0110012]. 

[36] E. Braaten and Y. Q. Chen, Phys. Rev. Lett. 76, 730 (1996) [arXiv:hep-ph/9508373]. 

[37] F. Yuan, C. F. Qiao and K. T. Chao, Phys. Rev. D 56, 1663 (1997) [arXiv:hep-ph/9701361]. 

[38] Y. J. Zhang, Y. Q. Ma, K. Wang and K. T. Chao, arXiv:0911.2166 [hep-ph]. 

[39] S. Fleming, A. K. Leibovich and T. Mehen, Phys. Rev. D 68, 094011 (2003) [arXiv:hep- 

ph/0306139]. 
[40] P. Pakhlov, arXiv:0901.2775 [hep-ex]. 

[41] Y. J. Zhang and K. T. Chao, Phys. Rev. Lett. 98, 092003 (2007) [arXiv:hep-ph/0611086]. 

[42] B. Gong and J. X. Wang, arXiv:0904.1103 [hep-ph]. 

[43] Y. Q. Ma, Y. J. Zhang and K. T. Chao, Phys. Rev. Lett. 102, 162002 (2009) [arXiv:0812.5106 
[hep-ph]]. 

[44] B. Gong and J. X. Wang, Phys. Rev. Lett. 102, 162003 (2009) [arXiv:0901.0117 [hep-ph]]. 

[45] Z. H. Lin and G. h. Zhu, Phys. Lett. B 597, 382 (2004) [arXiv:hep-ph/0406121]. 

[46] A. K. Leibovich and X. Liu, Phys. Rev. D 76, 034005 (2007) [arXiv:0705.3230 [hep-ph]]. 

[47] I. MaksA^myk, arXiv:hcp-ph/9710291. 

[48] E. Braaten and Y. Q. Chen, Phys. Rev. D 57, 4236 (1998) [Erratum-ibid. D 59, 079901 (1999)] 

[arXiv:hep-ph/9710357]. 
[49] C. B. Paranavitane, B. H. J. McKellar and J. P. Ma, Phys. Rev. D 61, 114502 (2000). 
[50] H. Jung, D. Krucker, C. Greub and D. Wyler, Z. Phys. C 60, 721 (1993). 
[51] W. L. Sang, L. F. Yang and Y. Q. Chen, Phys. Rev. D 80, 014013 (2009). 
[52] A. P. Martynenko, Phys. Rev. D 72, 074022 (2005) [arXiv:hep-ph/0506324]. 
[53] M. Gremm and A. Kapustin, Phys. Lett. B 407, 323 (1997) [arXiv:hep-ph/9701353]. 
[54] G. T. Bodwin, H. S. Chung, D. Kang, J. Lee and C. Yu, Phys. Rev. D 77, 094017 (2008) 

[arXiv:0710.0994 [hep-ph]]. 
[55] W. Y. Keung and I. J. Muzinich, Phys. Rev. D 27, 1518 (1983). 
[56] G. T. Bodwin and J. Lee, Phys. Rev. D 69, 054003 (2004) [arXiv:hep-ph/0308016]. 
[57] G. T. Bodwin and A. Petrelh, Phys. Rev. D 66, 094011 (2002) [arXiv:hep-ph/0205210]. 
[58] Y. Fan, Y. Q. Ma and K. T. Chao, Phys. Rev. D 79, 114009 (2009) [arXiv:0904.4025 [hep-ph]]. 
[59] Z. G. He, Y. Fan and K. T. Chao, arXiv:0910.3636 [hep-ph]. 

[60] E. Braaten and Y. Q. Chen, Phys. Rev. D 54, 3216 (1996) [arXiv:hep-ph/9604237]. 



41 



